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ABSTRACT 

At low temperatures, the main coolant in primordial gas is molecular hydrogen, H2. 
Recent work has shown that primordial gas that is not collapsing gravitationally but 
is cooling from an initially ionized state forms hydrogen deuteride, HD, in sufficient 
amounts to cool the gas to the temperature of the cosmic microwave background. This 
extra cooling can reduce the characteristic mass for gravitational fragmentation and 
may cause a shift in the characteristic masses of population III stars. Motivated by 
the importance of the atomic and molecular data for the cosmological question, we 
assess several chemical and radiative processes that have hitherto been neglected: the 
sensitivity of the low temperature H2 cooling rate to the ratio of ortho-H2 to para-fi2, 
the uncertainty in the low temperature cooling rate of H2 excited by collisions with 
atomic hydrogen, the effects of cooling from H2 excited by collisions with protons and 
electrons, and the large uncertainties in the rates of several of the reactions responsible 
for determining the H2 fraction in the gas. 

It is shown that the most important of neglected processes is the excitation of 
H2 by collisions with protons and electrons. This cools the gas more rapidly at early 
times, and so it forms less H2 and HD at late times. This fact, as well as several of 
the chemical uncertainties presented here, significantly affects the thermal evolution of 
the gas. We anticipate that this may lead to clear differences in future detailed three 
dimensional studies of first structure formation. In such calculations it has previously 
been shown that the details of the timing between cooling and merger events decides 
between immediate runaway gravitational collapse and a slower collapse delayed by 
turbulent heating. 

Finally, we show that although the thermal evolution of the gas is in principle 
sensitive to the ortho-para ratio, in practice the standard assumption of a 3:1 ratio 
produces results that are almost indistinguishable from those produced by a more 
detailed treatment. 



1 INTRODUCTION 



The very first stars to form in the Universe are be- 
lieved to have formed within small protogalactic objects, 
cooled primarily by molecular hydrogen (H2). Molecular hy- 
drogen cooling becomes ineffective at temperatures below 
T < 200 K, and at gas number densities n > 10 4 cm -3 , 
and so any dense fragments that form in the cooling 
and collapsing gas have a characteristic mass of a few 
hundred solar masses, set by the Jeans mass at this 
temperature and density ijAbel, Bryan, fc Normanl 120021 : 
Bromm, Coppi, fc Larsonl l2002j; see a lso the reviews of 
Bromm fc Larsonll2004l and lGloverll2005l k Since there is little 
evidence for sub-fragmentation during later sta ges of the co l- 
lapse (although for a dissenting view see lClark et al.l l2008). 
and since the high gas temperature leads to a high protostel- 
lar accretion rate, there seems little to limit the growth of 
the first stars, which may easily grow to ~ 100 Mr or more 
(see e.g. lYoshida et ailbOOrj ; lO'Shea fc Normanl |2007| ). 



Efficient cooling from hydrogen deuteride, HD, can al- 
ter this scenario. HD can cool the gas to lower tempera- 
tures than H2 , and remains an effective coolant up to higher 
densities, n ~ 10 6 cm -3 . The characteristic mass of stars 
formed in HD-cooled gas is therefore believed to be sma ller, 
~ 10 M (| Johnson fc Bromnj|2006l ; lYoshida et al.ll2007ft . re- 
flecting the smaller characteristic mass scale imprinted on 
the cooling gas. However, HD cooling will only bring about a 
change of this kind in the charact eristic mass scale if enough 
forms to cool the gas efficiently. iBromm. Coppi. fc Larsonl 
(2002) show that in simulations following the formation of 
the very first stars, in protogalaxies with virial tempera- 
tures T v i r < 10 4 K, this does not occur: the inclusion of 
deuterium chemistry and HD cooling has very little effect 
on the outcome. On the other hand, various authors have 
shown that in gas cooling from an initially ionized state, 
enough HD forms to cool the gas down to temperatures 
close to the temperature of the cosmic microwave back- 
ground (|Nakamura fc Umemurall2002l : iNagakura fc Omukail 
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2005; Johnson & Bromml 


2006;IShchekinov & Vasilievll2006; 
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20071; Yoshida et al J 12007). Note. 



however, that even without HD cooling the characteristic 
masses of objects collapsing from gas within a relic pri- 
mordial Hn region have al ready been demonstrated to be 
smaller (jO'Shea et al.ll2005h . 

This difference in thermal evolution, depending on 
whether or not the gas was once ionized, is a consequence of 
the chemistry of HD formation and destruction. The domi- 
nant reactions regulating the amount of HD in the gas are 



H 2 + D + 
and 

HD + H 4 



HD + H^ 



H 2 +D H 



(1) 



(2) 



Reaction 1 is exothermic, while reaction 2 is endothermic 
by 0.0398 eV (462 K), and so at low temperatures, chemi- 
cal fractionation occurs: the HD:H2 ratio becomes enhanced 
over the cosmological D:H ratio by a large numerical factor. 
Consequently, even though the HD cooling rate per molecule 
decreases with decreasing temperature, the HD cooling rate 
per unit volume can actually increase, owing to the in- 
crease in the HP abundance p roduced by this fractionation 
process (see e.g. iGlover! 120071 ). In conventional population 
HI st ar formation calculations (e.g. lAbel. Bryan, fc Normanl 
2002), the fractional ionization is small, and because of pdV 
heating the gas temperature never becomes low enough for 
chemical fractionation to become efficient. Therefore, HD 
cooling remains unimportant. In contrast, in gas cooling 
from an initially ionized state, more H 2 forms, owing to 
the non-equilibrium fra ctional ionization in the cooling gas 
jShapiro fc Kandll987h . and the gas can reach a lower tem- 
perature. In practice, the extra cooling provided by the en- 
hanced H 2 abundance is sufficient to cool the gas to a point 
at which chemical fractionation becomes very important, fol- 
lowing which HD dominates the cooling. 

Several processes and rate uncertainties, hitherto ne- 
glected, may interfere with this simple picture. First, most 
calculations assume a ratio of ortho-hydrogen (H 2 with nu- 
clear spin quantum number / = 1) to para-hydrogen (H 2 
with 7 = 0) that is (2/ ort ho + l)/(2/ par a + 1) = 3. This value 
is appropriate for warm H 2 in local thermodynamic equi- 
librium (LTE), which has many different rotational and vi- 
brational levels populated, but at low temperatures and low 
densities, the ortho-para ratio may differ significantly from 
this value. For instance, if only the J = and J — 1 rota- 
tional levels of the vibrational ground state are populated, 
then the equilibrium ortho-para ratio is 9 exp ( — 170.5/T). 
The relevance of this to the current situation lies in the 
fact that the energy associated with the v = 0, J = 2 ^ 
rotational transition in para-hydrogen, E20 = 509.85 K, is 
significantly smaller than the energy associated with the v = 
0, J = 3 — > 1 transition in ortho-hydrogen, E31 — 844.65 K. 
Consequently, para-hydrogen can cool the gas to lower tem- 
peratures than ortho-hydrogen. It is therefore possible that 
the ability of the gas to cool to the low temperatures re- 
quired for HD cooling to take over and dominate will be 
sensitive to the assumed ortho-para ratio, and that the out- 
come of calculations that determine it accurately will differ 
from that of calculations that assume a ratio of 3:1. 

A second issue affecting existing calculations is the 
fact that the low temperature behaviour of the H 2 cool- 



ing rate remains uncertain. The root cause of this uncer- 
tainty is the sensitivity of the low energy H-H 2 excitation 
cross-sections to the choice of the interaction potential used 
to calculate them. Most previous studies of HD cooling in 
primordial gas have made use of the fit t o the low-density 
H 2 cooling rate given by iGalli fc PalTa! il 19981). At T < 
600 K , this fit is based on excitation rates from Forrev et al.l 
1 19971 ) that wer e calculated using the B KMP2 potential en- 
ergy surface of |Boothrovd et al.l l|l996l ). However, recently 
IWrathmall fc Flower ( 2007 ) have published a new set of H 2 
collisional excitatio n rate coefficients based on calcula tions 
performed using the lMielke. Garrett, fc Pete rson (2002]) po- 
tential energy surface. The H 2 cooling function derived from 
these revised excita tion rates differs significantly from the 
IGalli fc Pailal i|l998h rate at temperatures T < 1000 K, but 
the consequences of this reduction in the cooling rate have 
yet to be explored in much detail. 

A third issue regarding the H 2 cooling rate is that fact 
that most previous calculations have only included the ef- 
fects of collisional excitation of H 2 by atomic hydrogen. How- 
ever, H 2 can also be excited by collisions with H 2 , He, H + 
and e~. As we show in t|3.2l in the conditions of interest 
for HD formation, several of these neglected processes play 
important roles. 

The final issue affecting studies of the role of HD cooling 
that we examine here is the impact of the large uncertain- 
ties that exist in several key rate coefficients for chemical 
reactions involved in the formation and destruction of H 2 . 
Although some of these uncertainties (which are discussed 
in detail in EI2.1I) ha ve re c eived previous study in the lit - 
craturc |Savin et all |2004| ; lei over. Savin fc J appscn 2006), 
their impact on the ability of the gas to cool to temperatures 
at which HD cooling becomes dominant has not previously 
been explored. 

In this paper, we explore these issues with the aid of 
a detailed chemical and thermal model of primordial gas, 
coupled to two simple dynamical models. Our main aim is 
to determine whether any of these sources of uncertainty 
can plausibly lead to significant differences in the evolution 
of the gas, or whether existing results on the role of HD 
cooling are robust. The structure of this paper is as follows. 
In Section [2l we outline the numerical model used in this 
work. In this context, we also discuss in more detail the ma- 
jor uncertainties highlighted above. In Section [3] we present 
and discuss our results, and we conclude in iJHwitli a brief 
summary. 



2 NUMERICAL MODEL 
2.1 Chemical network 

To model the chemistry of H 2 and HD in primordial gas, 
we use a chemical network consisting of 115 reactions be- 
tween 16 species, as summarized in Table |A"T1 This network 
differs significantly from previous treatments of primordial 
deuterium chemistry in that it includes the formation and 
destruction of doubly-deuterated hydrogen, D 2 . This is in- 
cluded because it has been suggested (D. Savin, private com- 
munication) that conversion of HD to D 2 at low gas temper- 
atures may be a significant destruction mechanism for HD, 
although in practice we find that it is unimportant. 
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For simplicity, we omit H^~, HeH + and their deuter- 
ated analogues from our chemical model. The abundances 
of these species are very small and their influence on the 
coolin g of the gas at intermediat e to low densities is mini- 
mal l|Glover fc S avin 20Q6|, I2008T ). so their omission should 
not significantly affect our results. We also omit lithium, for 
similar reasons. 

We assume that any radiation backgrounds are negli- 
gible and so do not include any processes involving pho- 
toionization or photodissociation. We also neglect cosmic 
ray ionization; the influence of this latter process on pro- 
moting HD cooling has bee n treated in detail elsewhere 
JShchekinov fc Vasilievl |2004 IVasiliev fc Shchekinovl 120061 ; 
IJasche. Ciardi fc Ensslinll2007l ; 1stacv fc Brommll2007l) . 

Whenever possible, rates for deuterated analogues 
of the basic hydrogen reactions have been taken from 
the primary literature, or fr o m the compilati o ns o f 
IStancil. Lepp fc Dakarnol (Il998h. I Wang fc Standi (|2002h 



and 
ever 



Walmslev. Flower fc Pineau des Foretsl l|2004 ). How- 



some reactions do not appear to have been previ- 
ously considered in the astrochemical literature. In cases 
where we have been unable to find an appropriate rate, 
we have generally adop t ed th e same procedure as in 
IStancil. Lepp fc Dalgarnol (|l998T l: for a non-deuterated re- 
action with a reaction rate that has a power-law tempera- 
ture dependence k oc T m , we have generated the rates of the 
deuterated analogues by multiplying this rate by a scaling 
factor (/xh/^d)" 1 , where fin and fio are the reduced masses 
of the reactants in the non-deuterated and deuterated reac- 
tions respectively. 

For reactions where the presence of a deuteron increases 
the number of distinguishable outcomes - e.g. the dissocia- 
tive attachment of HD with e~ (reactions 57-58), which can 
produce either H and D~ or H~ and D, in contrast to the dis- 
sociative attachment of H2 with e~ (reaction 23) which can 
only produce H~ and H - and where no good information 
exists on the branching ratio of the reaction, we assume that 
the probability of each outcome is uniform. For this particu- 
lar example, this gives branching ratios of 50% for reactions 
57 and 58 respectively. 

Finally, the rate coefficients for several of the included 
reactions require more detailed discussion, which can be 
found in sect ions 1 2 . 1.1 If2~l . 7 1 below . 

2.1.1 Associative detachment and mutual neutralization 
of H~ and D~ 

The rates of reactions 2 & 5, i.e. the associative detachment 
of H- with H: 

H^+H^Ha+e", (3) 
and the mutual neutralization of H~ with H + , 



H~ + H 4 



H + H, 



(4) 



are uncertain by up to an order of magnitude. When the 
fractional ionization of the gas is small, these uncertainties 
are unimportant, as in this case reaction 2 proceeds much 
faster than reaction 5. However, in gas with a high fractional 
ionization, such as gas recombining from an initially ionized 
state, reaction 5 competes with reaction 2 for the available 
H~ ions and so the uncertainties in the rates of these reac- 
tions introduce a significant uncertainty into the amount of 



H2 that is formed. A large associative detachment rate and 
small mutual neutralization rate lead to the production of 
a larger H2 fraction (at a given time) than a small associa- 
tive detachment rate and large mutual neutralization rate 
( Glo ver. Savin fc Jappsen|[2006l ). 

The default value for k 2 in our models is 



k 2 = 1.3 x 10~ 9 cm 3 s"\ 



(•») 



based on the measurement of ISchmeltekopf et al.l (| 1967t ) . 
However, in i|3.3l we present results from models performed 



(6) 



using 

k 2 = 5.0 x 10 -9 cm 3 s~ 
and 

k 2 = 0.65 x 10~ 9 cm 3 s 



(7) 



which repre sent plausible upper and lower bounds on the 
actual rate |Glover, Savin fc Jappsedl2006l ). 

Similarly, our default value for k$ is given by 



fc 5 = 2.4 x 10" 6 T 



-0.5 / 



1.0 + 



T 



20000 

taken from ICroft et all (|l999h . but in 
models using 



3 -1 
cm s 



(8) 

we also examine 



fc 5 = 5.7 x 10~ 6 T' - 5 + 6.3 x 10 -8 - 9.2 x 10 _11 T 5 

+ 4.4 x 10~ 13 Tcm 3 s~\ (9) 

taken from iMoselev et all (|l970l ) and 

fc 5 = 7.0xlO~V a5 cm 3 s" 1 , (10) 



taken 



from iDalgarno fc Leppl (|l987T ). 



I Glover. Savin fc Jappsenl (|2006l ) have suggested that 



the last of these rates may b e erroneously sm all, 



ing to typographical errors in IDalgarno fc Leppl (|l987T ) . 
Nevertheless, this rate has been used in a number of 
rece nt models of HD formati o n in primordial gas (se e 
e.g. iNagakura fc Qmukail 120051 ; Ijohnson fc Bromml 120061 ). 
justifying its consideration here. 

In view of the large uncertainties in the rates of reac- 
tions 2 & 5, we have assumed that identical rates apply for 
the deuterated analogues of these reactions (nos. 54-56 and 
66-68), since any small differences in the basic rates caused 
by the presence of one or two deuterons in place of protons 
are likely swamped by this basic uncertainty. 

2.1.2 Charge transfer from H + to H 2 (reaction 7) 

The most accurate cross-section for this process at astro- 
physic ally relevant energies is that computed by iKrstid 
l|2002t); the co r respon ding thermal rate coefficient is given 
in ISavin et~aH (|2004f ). However, as ISavin et~aH (|2004T ) dis- 
cuss in some detail, a large number of other rates for this 
reaction are given in the literature, differing by orders of 
magnitude at temperatures below 10 4 K. As this reaction 
is an important H2 destruction mechanism, particularly in 
gas recombining from an initially ionized state, and as most 
previous studies of HD formation in primordial gas have 
used one or another of these less accurate rate coefficients 
IIYoshida et al ] |2007i are a notable exception), it seems ap- 
propriate to examine the effect that the choice of this rate 
coefficient has on the final amount of H2 formed and on 



4 S. C. 0. Glover & T. Abel 



the ability of the gas to cool to temperatures at which 
HP cooling becom es dominant. Therefore, while we use the 
ISavin et al,l l|2004h rate in most of our models, we examine 
in i|3.3l the effect of using tw o other rates from the lit erature. 
The first of these, from lShapiro fe Kand (|l987T ) 



k 7 = 2.4 x l(T J exp | 



21200 ^ 3 _! 
— — J cm s 



(11) 



is, strictly speaking, onl y applicable to vibra t ionall y excited 
H2 , but in spite of this Jo hnson fc Bromml (|2006l ) use this 
rate for charge transfer with ground-state H2 in their study 
of HP cooling. As the comparison in Figure 1 of S avin et al.l 
(|2004T ) demonstrates, this rate is significantly larger than 

other determinations in th e literature . 

At the other extreme, lAbel et all fll997h quote a rate 

k 7 = exp (-24.2491469 
+ 3.4008244(ln T e ) 

- 3.8980040(lnT o ) 2 
+ 2.0455878(lnT c ) 3 

- 5.4161829 x lCr^lnTo) 4 
+ 8.4107750 x 10~ 2 (lnT c ) 5 

- 7.8790262 x 10~ 3 (lnT c ) 6 
+ 4.1383984 x 10" 4 (lnT c ) 7 

- 9.3634588 x 10" 6 (lnT c ) 8 ) cm 3 s" 1 , (12) 

where T c is the gas te mperature in units o f electron-volts. 
This rate is based on Ijanev et all (|l987f ), and has sub- 
sequently been adopted b y a number of authors (see e.g. 
Naga kura fc Om ukai 200i|). However, it is much smaller at 
T < 10 4 K than any of the other determinations in the 
ISavin et al.l (|2004l ) comparison. 



2.1.3 Collisional dissociation of H2 (reactions 8-11) 

In Table IA11 we list two rates for each process: one for H2 
that is all in the vibrational ground-state (appropriate for 
low density gas), and one for H2 with local thermodynamic 
equilibrium (LTE) level populations. At intermediate densi- 
ties, we adopt a rate coefficient for each reaction given by 



log h 



n/n c 



l + n/n c 



log h,_ 



1 + n/r 



logfc, v =o, (13) 



where ki is the collisional dissociation rate for collisions with 
species i, fc v =o,i and fci/rE,i are the rates for this reaction in 
the v = and LTE limits respectively, and n CI is the critical 
density, given by 



(14) 

n cr n C r,H n C r,H 2 ricr.He 

Here, in = riu/n, xu 2 ~ 2nn 2 /n, ihc = JiHc/", n is the 
number density of hydrogen nuclei, and 

n cr , H = dex [3.0 - 0.416 log T 4 - 0.327 (log T4) 2 ] , (15) 

n C r,H 2 = dex [4.845 - 1.3 log T 4 + 1.62 (logT 4 ) 2 ] , (16) 

and 

ncr.Ho = dex [5.0792 {l.O - 1.23 x 10~ 5 (T - 2000)}] , (17) 

with T 4 = T/IOO O O K. The expression for n cr , H is 
from iLepp fc ShulU (|l983l ). but has been decreased 



by an order of magn it ude, as recommended by 
iMartin. Schwarz fc Mandvl (Il996l) . The expression for 
ro C r,H 2 comes from IShapiro fc Kand (Il987l). and the expres- 
sion for n cri Hc comes from IPove et al.l (|l987l ). Note that 
this expression for the critical density assumes that in high 
density gas, n c <C riH, so that electron excitation of H2 does 
not significantly affect the value of n cr . 



2.1.4 He + recombination (reaction 19) 

In optically thick gas that is a mixture of neutral H and He, 
the effective He + recombination coefficient is given by 

fcig = 0.68fcl 9 , rr ,A + 0.32fci9, rr , B + fcig idi , (18) 

where we have assumed that nu 2 <S nu; see lOsterbrockl 
(1 19891 ) for a more detailed discussion. 

In these conditions, it is also necessary to take account 
of the photoionization of H caused by the He + recombination 
emission. As long as the gas is highly optically thick above 
the Lyman limit, this can be modelled as a local H ionization 
rate with a value 



Rpi = fc p in c n Hc+ cm 3 s , 
where 

fcpi = [0.68fcl9, rr ,A + 0.28fcl9,rr,B + fcl9,di] • 



(19) 



(20) 

We have not included a similar contribution from He ++ re- 
combination, as in conditions where the He ++ abundance is 
significant, we expect H to be almost completely ionized. 



2.1.5 Three-body H2 formation (reactions 30 & 31) 

At high densities, reactions 30 & 31 are important sources 
of H2. However, the rate coefficients for the se react i ons ar e 
highly uncertain, as previously discussed in iGloverl (|2007l ). 
To assess the importance of this uncertainty on our results, 
we have carried out simula tions using two different values 
for &3o: the first, taken from Abel. Bryan, fc Normanl (|2002l ) 
and partially based on lOrell ( 19871 ) is the lowest of the values 
we have found in the literature: 



1.14 x io- 31 T-°- 38 cm 3 s _1 
3.9 x 10- 30 T- 10 cm 3 s _1 



T< 300K 

T > 300 K 1 ' 



The o ther, taken from the recent paper of lFlower fc Harris! 
(2007) has the highest value at low temperatures of any of 
the rates we have found: 



fc 30 = 1.44 x io- 26 T -1 - 54 cm 3 s" 1 . 



(22) 



To fix the rate o f react ion 31, we follow 
iPalla. Salpeter. fc Stahleil (|l983h and assume that 

&31 = &30/8. 



2.1.6 Destruction of D2 by collision with H (reaction 107) 

The data tabulated in iMielke et all (|2003h span the temper- 
ature range 200 ^ T ^ 2200 K. At lower temperatures, we 
simply extrapolate our fit to the higher temperature data: 
this fit remains well behaved at low temperatures, and since 
the rate of this reaction falls off exponentially at low T, we 
are not particularly sensitive to errors in its value in this 
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temperature range. At T > 2200 K, we u se the simple expo- 
nential fit given by iMielke etHl (|2003l ) to their high tem- 
perature calculations; although not formally valid at these 
temperatures, the fit remains well-behaved, and hopefully 
lies not too far from the true value. 



2.1.7 Collisional dissociation of HD and D2 (reactions 
108-115) 

Fo r collisions with electrons, a ccura te rates are available 
in iTrevisan fc Tennyson! (|2002aT ) and iTrevisan fc Tennyson! 
i|2002bl ). For collisions with H, H2 or He, however, we have 
been unable to find a treatment in the literature. We have 
therefore assumed that the rates of these reactions in the 
v — and LTE limits are the same as for the correspond- 
ing H reactions (nos. 8-10). For D2, we also adopt the same 
value for the critical density, while for HD, we increase n cr 
by a factor of 100 to account for its larger radiative transi- 
tion probabilities. Note that although these rates are highly 
approximate, this probably does not introduce much un- 
certainty into the chemical model, as reactions 40 and 107 
become effective at much lower temperatures and therefore 
will generally dominate the destruction of HD and D2 in 
warm gas. 

2.2 The ortho-para hydrogen ratio 

In order to follow the evolution of the ortho-para hydro- 
gen ratio in the gas, we directly follow the time-dependent 
level populations of the lowest four energy levels of the H2 
molecule, the J = 0,1,2 and 3 rotational levels of the 
vibrational ground state. Rates for collisional transitions 
between these four states are taken from several sources: 
non-reactive collisions with H (which cannot change the 
ortho-para ratio) are treate d using the rates computed by 
IWrathmall fc Flowed i|2007l) , while for reactive collisions 
(which can c hange the ortho-para ra tio), we use the rates 
suggested by iLe Bourlot et all l|l999t) . Colli sions with pro - 
tons are treated using the rates computed by Gcrlich ( 1990h. 
Radia tive transitions rates are taken from IWolniewicz et al.l 
<| 19981 ). 

Newly-formed H2 is assumed, for simplicity, to reside in 
the J — ground state. This assumption is not correct: H2 
formed by associative detachment of H~ is, in general, highly 
excit ed and has a non-zero ortho-para ratio 1 Launav et al.l 
1991). However, it is easy to show that this assumption has 
little effect on the ortho-para ratio. In conditions where as- 
sociative detachment dominates the destruction of H~, we 
can write the H2 formation timescale as 

(23) 



10 _10 cm 3 s _1 , while from Table lATI we see that at a represen- 
tative low temperature of 200 K, ki — 2.07 x 10" 
Thus, the timescales are comparable only if 



^form — , 

K\X a - n 

where xh 2 an d x e - are the fractional abundances of H2 
and free electrons respectively. In comparison, collisions with 
protons cause the ortho-para ratio to reach equilibrium on 
a timescale 

fop 



too — 



(24) 



where f op is the ortho-para ratio, x K + is the fractional abun- 
dance of protons, and k op is an appropriately averaged rate 
coefficient for the co nversion of orth o-H2 to para-H2 by pro- 
ton collision. From iGerlichl l|l990h . we know that fc oc ~ 



xh 2 ~ 10 6 /op— 2 — 

X-H+ 



(25) 



If we make the reasonable assumption that x c ~ ~ x H + , and 
that fop is of order unity, then this argument demonstrates 
that the H2 formation process has a significant effect on the 
ortho-para ratio only when the H2 fraction is very small, 
^h 2 < 10~ 6 . 

A similar comparison can also be performed between 
the H2 formation timescale and the lifetimes of excited states 
of H2, but again the after-effects of the formation process are 
important only when the H2 fraction is very small. 

Our model for the ortho-para ratio becomes increas- 
ingly inaccurate at high temperatures, as the excitation of 
states with J > 3 or v > becomes important, but since 
the sensitivity of the H2 cooling rate to the ortho-para ra- 
tio is large only at low temperatures (see £12.31 below), this 
simplified approach is sufficient for our purposes. 



2.3 Thermal processes 

2.3.1 H2 cooling: collisions with H 

As we have already discussed in Section 1, the low tempera- 
ture rates for the collisional excitation of H2 by H are highly 
sensitive to the choice of potential energy surf ace used to de- 
scribe the H3 system (|Sun fc DalgarncTl 19941 ). An accurate 
determination of the H2 cooling function at low tempera- 
tures and low gas densities requires a level of accuracy in 
the potential that has been difficult to achieve, and as a 
consequence there are a number of determinations of the 
low-density limit of the H2 cooling function in the literature 
that differ substantially at temperatures T < 1000 K (see, 
for in stance, the comparison in Figure Al of iGalli fc PaHal 
1998). In recent years , the m ost widely used version has been 
that of lGalli fc Paiial (|l998h P1 

Ah 2 ,gp = dex [-103.0 + 97.59 log T - 48.05(log T) 2 

+ 10.80(log T) 3 - 0.9032(log T) 4 ] . (26) 

This rate is based on two separate sets of collisional rate 
coefficients. At tem peratures T < 600 K , the rates used are 
those computed bv lForrev et alj fl997) using a fully quan- 
tal approach and the BKMP2 potential energy surface of 
iBoothroyd et al.l (Il996l). A t T > 600 K, the rates used are 
those of lMandv fc Martini l|l993l ). which were computed us- 
ing the quasi-classical tra jectory approach and the LSTH 
potential energy su rface {Liul 11973 : ISiegbahn fcTiul Il978l : 
iTruhlar fc Horowitz|[l978l ). An ortho-para ratio of 3:1 is as- 
sumed at all temperatures. 

Recently, however, IWrathmall fc Flower! (|2007l ) 
have published a new set of collisional rate coeffi- 
cients computed using the pote ntial energy surface of 
IMielke. Garrett, fc Peterson! l|2002l ). The rms error in this 
new potential energy surface is more t han an order of 
magnitude smaller than the error in the IBoothroyd et all 



1 This cooling rate has units of erg cm 3 s 1 , as do all of the other 
cooling rates quoted in this paper, unless indicated otherwise. 
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Table 1. Fitting coefficients for the cooling rate of ortho-H2 ex- 
cited by collisions with atomic hydrogen 



Table 3. Fitting coefficients for the cooling rate of para-H2 ex- 
cited by collisions with H2 



Coefficient 


100 < T < 1000 K 


1000 < T < 6000 K 


Coefficient 


Para-H2 


Ortho-H 2 


Oo 


-24.330855 


-24.329086 


Oo 


-23.889798 


-23.748534 


01 


4.4404496 


4.6105087 


ai 


1.8550774 


1.76676480 


02 


-4.0460989 


-3.9505350 


0.2 


-0.55593388 


-0.58634325 


03 


-1.1390725 


12.363818 


03 


0.28429361 


0.31074159 


«4 


9.8094223 


-32.403165 


04 


-0.20581113 


-0.17455629 


or, 


8.6273872 


48.853562 


05 


0.13112378 


0.18530758 


06 


0.0 


-38.542008 








a 7 


0.0 


12.066770 









Table 2. Fitting coefficients for the cooling rate of para-H2 ex- 
cited by collisions with atomic hydrogen 



Coefficient 


100 < T < 1000 K 


1000 < T < 6000 K 


oo 


-24.216387 


-24.216387 


<n 


3.3237480 


4.2046488 


02 


-11.642384 


-1.3155285 


03 


-35.553366 


-1.6552763 


"4 


-35.105689 


4.1780102 


05 


-10.922078 


-0.56949697 


Oe 


0.0 


-3.3824407 


a 7 


0.0 


1.0904027 



(1996) potential, and IWrathmall fc Flowerl (|2007f) argue 
that it should allow a more accurate determination of 
the near-threshold behaviour of H2, and hence a bet- 
ter determination of the low-temperature excitation 



rates and cooling rate. Wrathmall fc Flowerl ([2007) and 



IWrathmall. Gusdorf fc Flowerl (|2007l ) show that there 



are indeed significant differences in the low temperature 
behaviour of a number of different excitation rates. 

We have used the r ate coefficients calculated by 
IWrathmall fc Flowerl (|2007h to compute separate cooling 
rates for ortho-Hh and para-H~2 in the low density limit due 
to collisions with atomic hydrogen. For ortho-H2, we find 
that the cooling rate in the temperature range 100 < T < 
6000 K is fit to within 2% with a function of the form 



7 

logA H2 ,H = CLi log(T3) 

i=0 



(27) 



where T3 = T/1000 K. The fitting coefficients a; are 
listed in Table [J Below 100 K, we extrapolate the 
IWrathmall fc Flowerl \200H ) rate as 

A oH2 ,h = 5.09 x 10- 27 T 3 1/2 exp (^^) . (28) 

Note that as HD cooling dominates at these low temper- 
atures, we are not particularly sensitive to errors in this 
extrapolation. 

For para-Fh, we follow a similar procedure: the para- 
H2 cooling rate for 100 < T < 6000 K can again be fit to 
within 3% by a function of the form of Equation 1271 using 
the fitting coefficients listed in Table [2] At T < 100 K, we 
use the extrapolation 

A pH2 ,h = 8.16 x W- 26 T^ exp (^^) . (29) 
Given these partial rates, the total H2 cooling rate in 



Table 4. Fitting coefficients for the cooling rate of ortho-H2 ex- 
cited by collisions with H2 



Coefficient 



Para-H2 



Ortho-H 2 



Oq 


-24.126177 


-24.020047 


ai 


2.3258217 


2.2687566 


02 


-1.0082491 


-1.0200304 


03 


0.54823768 


0.83561432 


04 


-0.33679759 


-0.40772247 


a 5 


0.20771406 


0.096025713 



the low-density limit due to collisions with atomic hydro- 
gen for gas with an ortho-hydrogen abundance x and para- 
hydrogen abundance x p is then simply 



At 



A H 2 ,H + 



A 



pH 2 ,H- 



(30) 



2.3.2 H2 cooling: collisions with H2 

To treat c ooling due to collisio ns between two H2 molecules, 
we follow Fl ower et all (|2000D and use rates for the excita- 
tion of pa ra-fi2 and ortho-f i 2 by ground-state para-Ii~2 de- 
rived from iFlower fc Rouefj | 1998) and rates for the excita- 
tion of par a-H~2 and ortho-Fb by g round-state ortho-H2 de- 
rived from lFlower fc Rouefil |l999). In the low density limit, 
and for temperatures in the range 100 < T < 6000 K, these 
rates are fit to within 2% by functions of the form 



logA H2 ,H 2 = 



5 

E 

i=0 



at log(T 3 y 



(31) 



The fitting coefficients are listed in Tables [3] and [4] for cooling 
from para-Ii~2 and ortho-H2 respectively. 

The total cooling rate in gas with an ortho-Hb abun- 
dance x and a para-Fb abundance x p is given by 



3^pApH 2 ,pH 2 "t" XpX ApH 2 ,oH 2 
"t" X Xp A H 2 ,pH 2 "f" ^qAoH 2 ,oH 2 



(32) 



where A p h 2 , p h 2 denotes the cooling rate due to the excita- 
tion of para-Hb by para-Ii~2, A p h 2 , h 2 the cooling rate due 
to the excitation of para-Hb by ortho-Hb, etc. 



2.3.3 H2 cooling: collisions with He 

The excitation of H2 by collisions with helium has been 
studied by a large number of authors (see e.g. iLee et all 
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Table 5. Fitting coefficients for the cooling rate of H2 excited by 
collisions with atomic helium 



Table 6. Fitting coefficients for the cooling rate of H2 excited by 
collisions with protons 



Coefficient 


Para-H2 


Ortho-H 2 


Coefficient 


Para-H2 


Ortho-H 2 


«o 


-23.489029 


-23.7749 


Oo 


-21.757160 


-21.706641 


ai 


1.8210825 


2.40654 


ai 


1.3998367 


1.3901283 


a 2 


-0.59110559 


-1.23449 


a 2 


-0.37209530 


-0.34993699 


13 


0.42280623 


0.739874 


0.3 


0.061554519 


0.075402398 


((4 


-0.30171138 


-0.258940 


04 


-0.37238286 


-0.23170723 


05 


0.12872839 


0.120573 


a 5 


0.23314157 


0.068938876 



120051 . and references therein). Currently, the most re- 
liable theoretical calculations appea r to b e those per- 
formed using the iMuchnick fc Russekl (Il994f) HeH 2 poten 
tial energy surface (e.g. IFlower. Roueff fc Zeippen 19981; 



1909 



Balakrishnan et al 



B alak rishnan. Forrev fc Dalgarnol. 

1999T) . The more recent iBoothrovd. Martin fc Peterson 
(2003) surface, which was expected to be more accurate, pro- 
duces results for some transitions t hat are in serious conflict 
with experimental determinations l|Lee et al.ll2005l ) and so 
results derived using this potential energy surface are cur- 
rently not considered reliable. 

In our models, we use an H2 cooling rate due 
to collisio ns with He that is deri v ed fr om the calcu- 
lations of IFlower. Roueff fc Zeippenl (|l998h for tempera- 
tures in the range 100 < T < 6 000 K and from 
Bala krishnan. Forrev fc Dalgarnol (ll999T ) for temperatures 
T < 100 K (which were not treated in the Flower et al . 
stu dy). Comparison of thelFlower. Roueff' fc Zeippenl (119981 ) 
and (Balakrishnan. Forrev fc Dalgarnol (|l999l ) rates at tem- 
peratures T > 100 K shows that they agree to within 10%. 
As before, we have derived separate rates for ortho-tb and 
para-Hb. In the low density limit, both cooling rates are fit 
to within 1% for temperatures T < 6000 K by a function of 
the form 



log A H2 ,He 



5 

£ 

i=0 



a, log(T 3 



(33) 



where T 3 = T/1000 K. The fitting coefficients en for the 
ortho and para cases are listed in Table [5] 

2.3.4 H2 cooling: collisions with protons and electrons 

In gas with a significant fractional ionization, collisions with 
protons and electrons can lead to a substantial H2 cooling 
rate. To treat the effect of collisi ons with pro tons, we use 
the rotational excitation r ates of i Gerlichl (|l990T ) and the vi- 
brational cross- sections of iKrstid l2002l ). In the low density 
limit, the cooling rates of ortho-Kb and para-H2 due to pure 
rotational transitions to levels with 2 ^ J ^ 7 in the vi- 
brational ground state, plus pure vibrational transitions to 
levels with 1 ^ v ^ 4 can be fit to within 2% over the 
temperature range 10 < T < 10000 K by a function of the 
form 



log A, 



iH+ = ^a 4 log(T 3 ) 

!=0 



(34) 



where T 3 = T/1000 K, using the fitting coefficients listed in 
Table [6] Note that these cooling rates include the effects of 
ortho-para interconversion in reactive collisions. In addition, 



A 



H 2 ,H+,0«1 



: 4.76 x 10" 



it is also necessary to account for the effect on the thermal 
balance of the gas of transitions from J = to J = 1 and vice 
versa. Conversion of para-H2 to ortho-H2 in the J = — ► 1 
transition cools the gas by 170. 5fc ~ 2.4 x 10~ 16 ergs per 
transition, while conversion of ortho-H2 to para-H2 in the 
J = 1 — > transition heats the gas by the same amount. In 
thermodynamic equilibrium, the number of transitions from 
J = to J = 1 exactly balances the number of transitions 
from J = 1 to J = 0, and so there is no net effect on the 
gas temperature. However, if the gas is not in thermody- 
namic equilibrium, then there can be net heating or cooling 
of the gas, depending upon whether the ortho-to-para ra- 
tio is greater than or less than the equilibrium value. We 
account for this in our model with a rate of the form 

9exp ( 1 t° 5 ) Xp ~ Xo ' ^ 

where x and x v are the fractional abundances of ortho-H2 
and para-H2 , and where we have again mad e use of the ro - 
tational excitation and de-excitation rates o f lGerlichl i|l990T l. 

We note that our treatment of H 2 cooling due to col- 
lisions with H + does not account for the effects of rov ibra- 
tional transitions, as the calculations by iKrstid (|2002T ) are 
not rotationally resolved. In view of the potential impor- 
tance of this process in primordial gas cooling from an ini- 
tially hot, ionized state, a more comprehensive treatment 
would be desirable. 

To treat H 2 excitation by collisions with free electrons, 
we us e the rates given by iDraine. Roberge fc Dalgarno 
□1983) , based on cross-sections from lEhrhardt et al. 
(1968), Cro mpton. Gibson fc Mcintosh! (I1969T) and 
iLinder fc Schmidt! 1 197lf ). IDraine. Roberge fc Dalgarnol 
( 19831 ) gives formulae for the collisional de-excitation 
rates of pure rotational transitions with AJ = 2 and 
pure vibrational transitions between v = 1,2 and 3 and 
the vibrational ground-state. Using these rates, we have 
computed the low density para-H2-e~ and ortho-H2-e~ 
cooling rates over a wide range of temperatures, and have 
fit them with functions 



1°S a h 2 , c - = log 



exp 



5 

x log(T 3 ) 



(36) 



where x = 509.85fc for para-H2 and x = 845fc for ortho- 
H 2 . The fitting coefficients are listed in Table [7] In both 
cases, the fit is accurate to within 10% over the tem- 
perature range 10 < T < 10000 K. We note that as 



IDraine. Roberge fc Dalgarnol l|l983l) do not gives rates for 
rotational transitions with A J > 2 or for rovibrational tran- 
sitions, our derived H2 cooling rates will underestimate the 
true rates at high temperatures. However, as collisions with 
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Table 7. Fitting coefficients for the cooling rate of H2 excited by 
collisions with electrons 



Coefficient 


Para-H2 


Para-H2 


Ortho-H 2 




T < 10 3 K 


T > 10 3 K 




do 


-22.817869 


-22.817869 


-21.703215 


01 


0.95653474 


0.66916141 


0.76059565 


02 


0.79283462 


7.1191428 


0.50644890 


03 


0.56811779 


-11.176835 


0.050371349 


14 


0.27895033 


7.0467275 


-0.10372467 


a 5 


0.056049813 


-1.6471816 


-0.035709409 



protons are considerably more effective at exciting Eh than 
collisions with electrons (see £12.3.51 below), the error that 
this introduces into our calculations is unlikely to be large. 

2.3.5 Hi cooling: the total cooling function 

In the low density limit, the total cooling rate per H2 
molecule (with units of erg s _1 ) is simply given by the sum 
of the cooling rates due to collisions with H, H2, He, H + and 
e~, i.e. 

A H2 ,n-»o = /J A H2 ,fcn k (37) 

k 

where k = H, H2, He, H + , e~. 

At high densities, the H2 level populations are in local 
thermodynamic equilibrium, and the H2 cooling rate per 
molecule is independent of the chemical composition of the 
gas and is given by 

Ah 2 ,lte = ^ AjiEjifj, (38) 

i,j>i 

where Aij is the radiative de-excitation rate for a transition 
from level j to level i, Eji is the corresponding energy, and fj 
is the fraction of H2 molecules in level j, co mputed assum- 
ing LT E. At intermediate densities, we follow iGalli fc Pallal 
l|l998l ) and write the H2 cooling rate as 



1 + AH 2 ,LTE/AH 2 ,n^0 

In Figure [T] we plot An 2 ,fc as a function of temperature 
for k = H,H2,He, H + and e~. In this plot, we assume a 
fixed ortho-para ratio of 3:1 (corresponding to x — 0.75 
and x p = 0.25). We also includ e in this plot, for the p urposes 
of comparison, the widely used lGalli fc Pallal (|l998f ) cooling 
function. 

It is immediately apparent from this plot that collisions 
between H and H2 are relatively ineffective at cooling the 
gas at low temperatures. Given equal abundances of H and 
H2, collisions with H2 provide more cooling than collisions 
with H for temperatures T < 1400 K. Similarly, collisions 
with He provide more cooling than collisions with H for T < 
2000 K, while collisions with protons or electrons are more 
effective over the whole of the temperature range examined 
here, again assuming equal abundances. 

Of course, in reality, the abundances of the various col- 
lision partners will generally not be equal: in low density pri- 
mordial gas, at the temperature of interest here, atomic hy- 
drogen is by far the most abundant species. Typically, in gas 
undergoing gravitational collapse within a small protogalaxy 




Temperature (K) 

Figure 1. H2 cooling rates per molecule, computed for n = 
10 — 4 cm -3 , for collisions with H (lower solid line), H2 (lower 
dashed line), He (dash-dotted line), e — (lower dotted line), and 
H+ (upper dotted line; note that this rate excludes the effects of 
transitions between the ortho and para ground states). In every 
case an ortho-para ratio of 3:1 is assumed. Transitions between 
the ortho and para ground states, brought about by collisions with 
H+, cool the gas at T> 150 K, and heat it at lower temperatures; 
note, however, that the low-temperature heating is a consequence 
of our adoption of a tempera ture-independent ort ho-para ratio. 
Also shown is the widely used Icialli & Palla (1998) cooling func- 
tion (upper solid line), which considers only collisions between H 
and H2. 



(Tvir < 10000 K) forming in a region of the intergalactic 
medium (IGM) not yet affected by stellar feedback, one finds 
abundances relative to atomic hydrogen of ihc = 0.0825 for 
He, xn 2 ~ 10 -3 for H2 and x u + ~ x e - ~ 10~ 4 for pro- 
tons and electrons. If these relative abundances are taken 
into account, then atomic hydrogen becomes comparatively 
more effective. Collisions with H2 become completely unim- 
portant for the whole of the temperature range studied, and 
collisions with electrons can also be neglected. However, col- 
lisions with helium remain important at low temperatures, 
and in fact dominate the H2 cooling rate for T < 650 K, 
despite the significantly larger abundance of hydrogen rel- 
ative to helium. Collisions with protons are also important 
at T < 400 K, in spite of the low proton abundance. 

In gas cooling from an initially ionized state, similar 
conclusions hold regarding the relative importance of col- 
lisions with H, H2 and He. However, in this case, values 
for x u + and x e - that are 10-100 times larger are not un- 
common, and the effects of H2-H + and H2-e~ collisions are 
therefore much greater. 

It is also interesting to compare the r elative impor- 
tance of the various processes if one adopts the lGalli fc Pallal 
(1998) rate for cooling from H2-H collisions in pl ace of 
our value derived from IWrathmall fc Flowerl (|2007). Fig- 
ure \T\ demonstrates that the lGalli fc Pallal cooling rate pro- 
vides significantly more co oling at T < 1000 K than the 
newer IWrathmall fc Flower! cooling rate, with the rates dif- 
fering mo st significantly at temperatures 300 < T < 500 K, 
where the lGalli fc Pallal rate provides almost five ti mes more 
cooling than the comparable IWrathmall fc Flowerl rate. Be- 
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cause of this, collisions with He and with protons and elec- 
trons are less effective in compar ison to collisio ns with H at 
T < 1000 K when one uses the iGalli fc Pallal rate. Never- 
theless, even though it is no longer the dominant process, 
cooling from Eb-He collisions remains important at low tem- 
peratures, as it can contribute 20-30% of the total H2 cool- 
ing rate. Furthermore, H2-H + and H2-e~ collisions will also 
still be important if the fractional ionization of the gas is 
large (:r H+ > 10" 3 ). 

We should note at this point that we are not the 
first authors to highlight the potential importance of 
H2-He and H2-H + col lisions for cooling primordial gas. 
iLe Bourlot et all (| 19991 ) include the effects of Hh-He col- 
lisions in their calcula t ions of the H2 cooling function, 
as do ISantoro fc Shulll (|2006t ): the i mportance of helium 
is also discussed at some length in iFlower et al.l (|2000t ). 
Th e possible im p ortan ce of H2-H + collisions was noted 
by IGalli fc Pallal (|l998l) and their effects were examined 
in more deta i l by Flower fc Pineau des Foretsl (|2000h and 
IFlower et al.l (|200C ). although only pure rotational transi- 
tions were considered. On the other hand, to the best of our 
knowledge, we are the first authors to consider the effects of 
H2-e~ collisions in primordial gas. 

Finally, although in this section we have given fits to 
the cooling rates of ortho-Eh and para-H2 separately, since 
we are interested in the effects of varying the ortho-para 
ratio, we recognize that for some purposes it may be useful 
to have the rates for a gas that has the often-assumed 3:1 
mix of ortho and para-Hb. In Table [8] we list fits to the low 
density H2 cooling rates due to collisions with H, H2, He, 
H + and e~ for this case. All of these fits are of the form 
5 

logA H2 ^^adogtTa) 8 , (40) 

i=0 

where T 3 = T/1000 K and the accuracies are comparable 
to the accuracies of the separate ortho and para-H2 fits. 
Note that the H2-H + rate quoted here does not include the 
effects of collisional transitions from J = to J = 1 or 
vice versa. However, this can be included through the use 
of Equation 1351 with x p — 0.25 and x — 0.75 for the 3:1 
ortho-para ratio case. 

2.3.6 H2 cooling: sensitivity to the ortho-para ratio 

In Figure [2] we compare three different H2 cooling rates: one 
for pure ortho-H2, one for pure para-H2, and one for which 
we assumed the standard 3:1 ortho-para ratio. In each case, 
we assume that n <C n C rit, so that we are in the low density 
limit, and adopt fractional abundances relative to hydrogen 
of iHe = 0.0825, xn 2 — 0.001, and x H + — x c - = 10 -4 for 
He, H2, H + and electrons respectively. Note that at T < 
230 K in the ortho-H 2 case and for T < 98 K in the 3:1 
ratio case, collisional conversion of ortho-H2 in the J = 1 
rotational level to para-H2 in the J = rotational level 
by protons heats the gas, and that the lowest temperature 
portions of the curves plotted in Figure[5]for these two cases 
therefore represent heating rates. 

The figure demonstrates the importance of the ortho- 
para H2 ratio in determining the H2 cooling rate at temper- 
atures below a few hundred K. For instance, at T — 300 K, 
there is a difference of an order of magnitude between the 



Table 8. Fitting coefficients for H2 cooling rates, for a 3:1 ortho- 
para ratio 



Species Temperature range (K) Coefficients 



H 


10 < T < 100 


a = -16.818342 
ai = 37.383713 
a 2 = 58.145166 
a 3 = 48.656103 
a 4 = 20.159831 
as = 3.8479ol0 


H 


100 < T ^ 1000 


a = -24.311209 
01 = 3.5692468 
a 2 = -11.332860 
a 3 = -27.850082 
a 4 = -21.328264 
£15 = — 4.2olyUzo 


H 


1000 < T ^ 6000 


a = -24.311209 
ax = 4.6450521 
a 2 = -3.7209846 
a 3 = 5.9369081 
a 4 = -5.5108047 
05 = 1.5538288 


H 2 


100 < T < 6000 


a = -23.962112 
01 = 2.09433740 
a 2 = -0.77151436 
a 3 = 0.43693353 
a 4 = -0.14913216 
as = — 0.033D38326 


He 


10 < T < 6000 


a = -23.689237 
ai = 2.1892372 
a 2 = -0.81520438 
a 3 = 0.29036281 
a 4 = -0.16596184 
05 = 0.19191375 


H+ 


10 < T < 10000 


a = -21.716699 
01 = 1.3865783 
a 2 = -0.37915285 
a 3 = 0.11453688 
a 4 = -0.23214154 
a 5 = 0.058538864 


c~ 


10 < T < 200 


a = -34.286155 
ai = -48.537163 
a 2 = -77.121176 
a 3 = -51.352459 
a 4 = -15.169160 
a 5 = -0.98120322 




200 < T < 10000 


a = -22.190316 
ai = 1.5728955 
a 2 = -0.21335100 
a 3 = 0.96149759 
a 4 = -0.91023195 
a s = 0.13749749 



cooling rate of para-H2 and the cooling rate of ortho-H2. 
Because of this large disparity, para-H2 will provide most 
of the contribution to the H2 cooling rate at these temper- 
atures even in gas that contains primarily ortho-H2; e.g. it 
is the 25% of para-H2 that provides most of the cooling in 
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Figure 2. Comparison of the H2 cooling rate per H2 molecule, 
for gas with only ortho-H2 (solid line), para-H2 (dashed line) or 
a 3:1 mix of ortho and para-H2 (dotted line). We assume that 
xhc = 0.0825, xn 2 = 0.001 and x H + = x e - = 10 -4 . Note that 
for T < 230 K in the ortho-H 2 case and T < 98 K in the 3:1 ratio 
case, the rate plotted is the net heating rate, after accounting for 
heating due to the collisional conversion of J = 1 ortho-H2 to 
J = para-H2 by protons. 



the 3:1 case. Consequently, relatively small deviations in the 
ortho-para ratio may have a large effect on the low temper- 
ature H2 cooling rate. 

Furthermore, although the specific values for the cooling 
rates plotted in Figure [2] are sensitive to our assumed chem- 
ical abundances, the basic point that the low temperature 
H2 cooling rate is highly sensitive to the assumed ortho- 
para ratio is robust, as it is a consequence of the difference 
in the energy separations of the lowest levels of para-H2 
(E 20 = 509.85 K) and ortho-H 2 (£31 = 844.65 K). 



2.3.7 HD cooling 

To model HP cooling, we use the coolin g function of 
iLipovka, Nufiez-Lopez. fc Avila-Reesd l|2005h . This param- 
eterization of the HD cooling rate assumes that HD-H col- 
lisions make the dominant contribution. This is a much 
safer assumption in the case of HD cooling than in the case 
of H2 cooling. Excitation rate coefficients for HD-H colli- 
sions are typicall y muc h larger than for H2-H collisions, 
and Hi ower et ah! (|2000l ) show that they are comparable 
to the excitation rates for HD-He or HD-H2 collisions. As 
«h 2> tihc S> nu 2 in the conditions of interest here, this 
means that the HD-H contribution will dominate. 

The larger excitation rates for HD-H collisions also re- 
duce the importance of collisions with protons or electrons. 
Although accurate excitation rate coefficients for HD-H + or 
HD-e~ collisions do not appear to be available, it seems rea- 
sonable to assume that they will be of a similar order of 
magnitude to the corresponding processes with H2. If so, 
then in the T < 200 K temperature regime in which HD 
cooling is important, collisions with electrons or protons be- 
come comparable to collisions with atomic hydrogen only 
for fractional ionizations x > 0.1. To find such a large frac- 
tional ionization in gas this cold would appear to be highly 



unlikely, and so it seems relatively safe to neglect the effects 
of collisions with pro tons or electrons. 

Although the ILipovka. Nunez-Lopez, fc Avila-R ccsc 

(2005) parameterization of the HD cooling rate is formally 
valid only in the temperature range 100 < T < 2 x 10 4 K, 
we have compared its behaviour at lower temperatures 
with an explicit calculation of the cooling rate made using 
radiat ive de-excitation rates from lAbgrall. Roueff fc Vialal 
( 1982) an d collisional rates extrapolated from those com- 
puted b y IWrathmall. Gusdorf fc Flowed (I2007T ) . We find 
that the ILipovka. Nunez- Lopez, fc Avila-ReeseT(l2005l ) rate 
remains reasonably accurate down to temperatures as 
low as 50 K, with errors no greater than 20%, and 
that even at T = 30 K it remains accurate to within 
a factor of two. At temperatures T 100 K, the 

ILipovka. Nunez-Lopez, fc Avila-Reesd (|2005h cooling rate 
slightly underestimates the e ffects of HD cooling compared 
to the newer calculations of Wr athmall. Gusdorf fc Flowerl 
|2007l ). presumably owing to the more accurate vibrational 
excitation rates used in the latter, but the differences are 
relatively small and in any case occur in the temperature 
reg ime in which H2 cooling dominates. The bre akdown of 
the lLipovka. Nunez-Lopez, fc Avila-Reesel ([2005) fit at very 
high temperatures (T > 20000 K) is unimportant, as the 
gas in our models never exceeds this temperature, nor could 
HD cooling ever be significant in this temperature regime 
where the cooling from the Lyman a line of neutral hydro- 
gen peaks. 

To correctly model the effects of HD cooling at low tem- 
peratures, it is necessary to take the effects of the cosmic 
microwave background (CMB) into account. We do this ap- 
proximately, by using a modified HD cooling rate, A H d, de- 
fined as 



Ahd = Ahd(T) — Ahd (Tomb) 



(41) 



where Ahd(T) and Ahd (Tomb) are the unmodified HD cool- 
ing rates at the gas temperature T and the CMB tempera- 
ture Tcmb respectively. 

The quoted range of densiti es for which the 

ILipovka. Nunez-Lopez, fc Avila-Reesd l|2005h cooling func- 
tion is valid is 1 < n < 10 s cm -3 . To extend the range of the 
cooling function to densities n < 1 cm -3 , we assume that 
at these densities the cooling rate per molecule is directly 
proportional to n, and hence that 

Ahd(?i — n) — nAuo{n = 1) (42) 

for n ^ 1 cm -3 , where Ahd(tc) is the cooling rate per HD 
molecule (with units ergs -1 ) at gas number density n. To 
extend the cooling function to high densities, n > 10 8 cm -3 , 
we assume that the HD molecule is in LTE and thus has a 
cooling rate per molecule that is independent of density. In 
this regime, 



A HD (n > 10 s ) = Ahd (n = 10 s ). 
In view of the fact that 1 <C n C r,HD <C 10 s 



(43) 



cm , where 

ficr.HD is the HD critical density, both of these assumptions 
appear well justified. 



2.3.8 HJ, HD+ and DJ cooling 

In view of its possi ble importance in hot, ionized gas (see 
lYoshida et al.l 120071 ) , we include the effects of vibrational 
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cooling from the molecular ion, as well as from its deuter- 
ated analogues HD + and D^. 

At low densities, the main contributions to the cool- 
ing rate come from excitatio ns by collisions with electron s 
and with neutral hydrogen (|Suchkov fc Shchekinovl 1 1978ft . 
We have computed the cooling rate d ue to collisions with 
electro ns, using the vibrational rates of ISarpal fc Tennyson! 
\ 19931 ) for excitations from v = to v = 1 and v = 2; exci- 
tations to higher vibrational states (v = 3 to v = 8) are also 
included, under the assumption that the de-excitation rates 
for these transitions are comparable to the de-excitation rate 
from v = 2 to v — 0. We note that even at temperatures as 
high as 10 4 K, at least half of the total cooling comes from ex- 
citations to v — 1 and v — 2, and so our approximate treat- 
ment of transitions to the higher vibrational states makes 
the cooling rate uncertain by at most a factor of a few at 
high temperatures (and by far less than this at low temper- 
atures). The resulting cooling rate, A cH +, is given at 
T sC 2000 K by 

3025 \ 



A 



1 1 1 n — 19rri — 0.34 

,+ = 1.1 x 10 1 exp I 



and at T > 2000 K by 



A „+ = 3.35 x 10" 21 T° 12 exp | 



T J 



3025 s 



T J 



(44) 



(45) 



n — c,H 



+ n c - + A, 



,+ tlH. 



(50) 



To model cooling from vibrational transitions in HD + , 
we assume, in the absence of better information, that the 
low density cooling rate is the same as that used for Hj. 
However, since HD + has much larger radiative transition 
rates than Hj, the LTE cooling rate for HD + is much larger 
than that for Hj. We have calculated the HD + LTE c ooling 



rate per ion using level ene rgies from Karr fc Hilicol (12007 ) 
and t ransition rates from IPeek. Hashemi- Attar, fc Beckell 
(1979), and have fit it with the function 



= 1.09 x 10 _11 T - 03 exp 



2750 



T J 



at temperatures T ^ 1000 K and 



A 



LTE.HD + 



5.07 x 10" 12 T°- 14 exp 



2750 



T J 



(51) 



(52) 



at T > 1000 K. For densities between the low density and 
LTE limits, we use a function of the form of Equation [35] to 
compute the HD + cooling rate. 

Finally, to model cooling, we simply assume that 
the same rates apply as for Hj cooling. In practice, the very 
small size of the typical D J abundance renders this process 
irrelevant. 



T s: 1000 K a rate 



A H H + = 1.36 x 10~ z " exp | 



and at T > 1000 K a rate 



3152 \ 
T ) 



(46) 



A 



h,h; 



= dex [-36.42 + 5.95 log(T) - 0.526 log(T) 2 ] . (47) 



The h igh temperature rate is a fit made by Galli fc Palla 
l|l998ft to the rate given in ISuchkov fc Shchekinovl i| 19781 ); 
note that owing to a normalization e rror, the rate s given 
for A„ „+ and A „+ in Figure A2 of iGalli fc PaUal (1 19981 ) 

rl,n 2 e ' M 2 

are too large by a factor of ten. The low temperature rate 
given here is a physically re asonable extrapolation of the 
ISuchkov fc Shchekinovl (jl978l ) rate that has the correct ex- 
ponential fall-off at low temperature. 

At high densities, the vibrational levels of Hj will be 
in LTE. In this regime, the cooling rate per Hj ion is given 
approximately by 



At 



,+ = 2.0 x 10" 



-19^0.1 

1 exp i 



3125 \ 
T ) 



(48) 



To compute this rate, we included contributions from all 
vibrational st a tes v ^ 8 and used level energies from 
Karr fc Hilicol (l2007f) and radiative transition rates from 



Posen. Dakarno fc Peek] (j 19831 ). The effects of rotational ex- 



citation were not included, but are unlikely to change this 
expression by a large amount, owing to the very small tran- 
sition rates associated with these transitions. 

At intermediate densities, we assume that the Hj vi- 
brational cooling rate per Hj ion is given approximately by 
the function 

A T 



A H+ 



LTE.H+ 



1 + A LTE,H+/ A n-,0,H+ ' 



(49) 



where A n _ tQ H + is the cooling rate per Hj ion in the low 
density limit, and is given by 



2.3.9 Other processes 

In addition to the coolants listed above, we also in- 
clude radiative cooling from the electro nic e x citati on of 
H, He and He + using rates taken from ICenl l|l992l ) and 
Bray et all (|2000l ). Compton cooling (again using a rate from 
Cenlll992l). and brem sstrahlung (using the rates given in 
Shapiro fc Kandl 1987ft . 

Moreover, we also include the effects of chemical cooling 
from the collisional ionization of H, He and He + (reactions 
12, 17 and 18), collisional dissociation of H2 (reactions 8- 
11), the destruction of H2 by charge transfer (reaction 7), 
and the recombination of H + , He + and He ++ (reactions 13, 
19 and 20), as well as chemical heating arising from the 
formation of H2 via reactions 2 and 4. 

Furthe r details of our treatment of these processes can 
be found in iGlover fc Jappser] (|2007j ). 



2.4 Model setup and initial conditions 

We model the chemical and thermal evolution of primordial 
gas within the context of two simple toy models for its dy- 
namical evolution. In one, we assume that the gas evolution 
is isobaric. In this model, the gas temperature evolves as 



dT 
At 



7 



- 1 (j, a - r 

7 k p 



r dln/x 
dt 



T dln7 

' - 1 dt 



(53) 



where [i is the mean molecular weight of the gas (in 
grammes), A and V are the total cooling and heating rates 
per unit volume, and the other symbols have their usual 
meanings. The rate of change of [i can be easily determined 
from the chemical rate equations. For the adiabatic index 7, 
we use the expression 

5 + 5xhb + 5x e — 3xn 2 



7 ; 



3 + 3XHe + 3x e - lH 2 



(54) 
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where ihc, xh 2 an d x c are the fractional abundances of he- 
lium, H2 and free electrons relative to the abundance of hy- 
drogen nuclei. In practice, xn 2 remains small at all densities 
encountered in our isobaric models, and so 7 ~ 5/3 through- 
out. Finally, at any point in the evolution of the gas, we can 
relate the gas density to the temperature by 



(?) 



/'i 



(55) 



where pi , TS and pi are the initial values of the density, tem- 
perature and mean molecular weight respectively. 

In the other model, we assume that the gas undergoes 
gravitational collapse at the free-fall rate. In this model, the 
gas temperature evolves as 



dT _ 7-1 
dt ~ p 



+ ^^ + TiiM (56) 



7 - 1 dt 



dt 



and the gas density evolves as 



dp 
dt 



£_ 

iff 



(57) 



where ta — -\/37r/32Gp is the free-fall time. 

The first of these models approximates the case of gas 
that has been strongly shocked (e.g. by a supernova blast- 
wave or in a halo merger) but that is not yet gravitation- 
ally unstable. The free-fall collapse model approximates the 
other extreme case, in which gas is highly gravitationally un- 
stable, and contracts at the maximal rate. Realistically, the 
dynamical evolution of gas involved in high-redshift struc- 
ture formation probably lies somewhere in between these 
two cases. 

Although one could use far more sophis ticated models 
for th e dynamical evolution of the gas (see e.g. lYoshida et ail 
2007), our use of these simple models allows us to rapidly 
explore the effects of the various different sources of uncer- 
tainty discussed in this paper, and to highlight which are 
deserving of more numerically expensive three dimensional 
studies, and which are unimportant and can be safely ig- 
nored in future work. 

To evolve the coupled set of chemical rate equations and 
the thermal energy equation we us e the DVODE solver of 
iBrown. Byrne, fc Hindmarshl (|l989h - 

We adopt standard helium and deuterium abundances 
of ZHe = 0.08 25 and id = 2.6 x 10~ 5 relative to hydrogen 
l|Molardl200"7l ). and begin our simulations with gas in which 
hydrogen and deuterium are fully ionized and the helium 
is singly ionized. The initial abundances of all other species 
are set to zero. 

We fix the initial temperature at T\ = 20000 K, and ex- 
amine models with three different initial densities: m = 0.03, 
1, and 30 cm -3 . We run all of our models for two different 
redshifts, z = 10 and z = 20; the latter value is perhaps 
more appropriate for the study of the earliest generations 
of star formation, but the former allows us to look at the 
effects of having a CMB temperature that is much smaller 
than the temperature that the gas can reach through H2 
cooling alone. 



3 RESULTS 

3.1 Ortho-para ratio 

In order to establish the effect that variations in the H2 
ortho-para ratio have on the thermal evolution of primor- 
dial gas, we considered four separate cases: the two limit- 
ing cases in which all of the H2 is in the form of ortho- or 
para-hydrogen respectively, a third case in which the stan- 
dard ratio of 3:1 was assumed, and a final case in which the 
ortho-para ratio was determined self-consistently, although 
approximately, from the populations of the lowest four ro- 
tational levels, as outlined in £12.21 

In Figure [3^,, we show how the gas temperature evolves 
as a function of density for these four cases in two free- 
fall collapse models with initial density n\ — 0.03 cm -3 for 
redshifts z = 10 (lower set of curves) and z — 20 (upper set 
of curves). Figures [3p and show similar results for models 
with initial densities m = 1 and 30 cm -3 respectively. In 
each figure, the dashed and dash-dotted curves correspond 
to calculations in which the H2 is all in ortho or para form 
respectively, the solid curves correspond to the calculations 
that assume an ortho-para ratio of 3:1 and the dotted curves 
correspond to the calculations in which the ortho-para ratio 
was determined dynamically. The horizontal dashed lines 
give the CMB temperature at z = 20 (upper line) and z = 10 
(lower line). 

Figure |3] demonstrates that there are significant differ- 
ences between the temperature evolution in the pure ortho- 
H2, pure para-H2, and 3:1 ratio calculations. The para-Eb 
and 3:1 ratio calculations differ primarily at n < 10 5 cm -3 , 
with temperatures differing by as much as 50% at n ~ 
100 cm -3 . At n > 10 5 cm -3 , however, the simulations be- 
come convergent, and little difference remains in the tem- 
perature evolution. The ortho-Hb simulations show an even 
greater difference in behaviour. In most of these simulations, 
the gas temperature remains significantly larger than in the 
para-Eb or 3:1 ratio runs, differing by a factor of two or more, 
and failing to converge with the other simulations even at 
n > 10 5 cm" 3 . The gas temperature in most of the ortho-Eb 
simulations remains above 100 K throughout the collapse, 
and although a comparison of cooling rates shows that HD 
cooling does become dominant in these simulations, it does 
not succeed in driving down the temperature to the same 
extent as in the other runs. The one exception is the sim- 
ulation with ni = 0.03 cm -3 and z = 10, which does cool 
significantly below 100 K and which converges with the cor- 
responding para-H2 and 3:1 ratio calculations. 

Despite the apparent sensitivity of the temperature evo- 
lution to the ortho-para ratio, Figure [3] demonstrates that 
there is essentially no difference between the results of cal- 
culations in which the ortho-para ratio is fixed at 3:1 or 
calculated self-consistently from the H2 level populations. 
Figure [4] helps to demonstrate why this is so. In the figure, 
we show the dependence of the ortho-para ratio on the gas 
temperature in a representative free-fall collapse model with 
z = 20 and m = 1 cm -3 . At temperatures T > 200 K, the 
ortho-para ratio is approximately three, both at low den- 
sities (solid line, lower branch) and at high densities (solid 
line, upper branch). At lower temperatures, the ortho-para 
ratio falls off steeply with decreasing temperature, and at 
the lowest temperature reached by the gas, T ~ 81 K, para- 
hydrogen is almost as abundant as ortho-hydrogen. How- 
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Figure 3. (a) Temperature evolution as a function of gas num- 
ber density in free-fall collapse models with initial density n\ = 
0.03 cm -3 and initial redshifts z = 20 (upper set of curves) 
and 2 = 10 (lower set of curves). Four cases are examined: gas 
with pure ortho-ti2 (dashed curves), pure para-H2 (dash-dotted 
curves), a 3:1 ortho-para ratio (solid curves) and an ortho-para 
ratio determined by solution of the simplified level population 
calculation discussed in £12.21 (dotted curves; note that these are 
barely distinguishable from the solid lines in the plot). The CMB 
temperature at z = 10 and z = 20 is indicated by the horizontal 
dashed lines, (b) As (a), but for gas with rij = 1cm -3 . (c) As (a), 
but for gas with n\ = 30 cm -3 . 
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Figure 4. Ortho-para ratio as a function of gas temperature in 
a free-fall collapse model with z = 20 and n; = 1 cm -3 . The 
collapsing gas evolves initially towards lower temperatures along 
the lower branch of the solid line, before progressing back toward 
higher temperatures along the upper branch as the gas heats up 
at high density. The dotted line gives the equilibrium ortho-para 
ratio for our simplified model H2 molecule. The vertical solid 
line indicates the temperature below which HD cooling becomes 
dominant. 



ever, in this calculation, HD cooling dominates over H2 cool- 
ing at a temperature T = 135 K (indicated in the figure by 
the vertical solid line). At this temperature the ortho-para 
ratio is ~ 2.4, and so the H2 cooling rate does not differ 
greatly from the rate that we obtain by assuming a fixed 
ortho-para ratio of 3:1. In other words, at the temperatures 
where the true H2 cooling rate differs significantly from the 
H2 cooling rate in the 3:1 ortho-para case, H2 cooling is itself 
unimportant, and HD cooling dominates. We have verified 
that the same explanation also serves to explain the results 
of our other free-fall collapse models. 

Finally, to check that our conclusions do not depend on 
our choice of dynamical model, we have examined the be- 
haviour of isobarically evolving gas in the same four cases, 
as illustrated in Figure[5] We see that again the temperature 
evolution is sensitive to extreme variations of the ortho-para 
ratio, but that the results of the self-consistent calculation 
are barely distinguishable from those of the calculation as- 
suming a 3:1 ortho-para ratio. Further investigation demon- 
strates that the reason for this similarity is the same as in the 
free-fall collapse case: at temperatures where H2 cooling is 
significant, the ortho-para ratio remains close to three, while 
at the low temperatures at which it differs significantly from 
three, H2 cooling is unimportant and HD cooling dominates. 

We can therefore conclude that the adoption of a fixed 
ortho-para ratio of 3:1, although strictly speaking unjustified 
at T < 200 K, is nevertheless an adequate assumption for 
modelling the temperature evolution of primordial gas, and 
that therefore this potential source of uncertainty ultimately 
proves to be unimportant. 
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Figure 5. (a) Temperature evolution as a function of time in a 
set of models in which the gas evolution is isobaric. These models 
have an initial gas number density n; = 0.03 cm -3 , and results 
are plotted for both z = 10 (lower set of lines) and z = 20 (upper 
set of lines). Four cases are examined: gas with pure ortho-H2 
(dashed lines), pure para-H2 (dash-dotted lines), a 3:1 ortho-para 
ratio (solid lines) and an ortho-para ratio determined from the 
level populations of our model H2 molecule (dotted lines; note 
that again these are barely distinguishable from the solid lines 
in the plot). The CMB temperature at z = 10 and 2 = 20 is 
indicated by the horizontal dashed lines, (b) As (a), but for an 
initial number density n\ = 1 cm -3 , (c) As (a), but for an initial 
number density rij = 30 cm~ 3 . 
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Figure 6. Temperature evolution as a function of gas num- 
ber density in free-fall collapse models with initial density n; = 
1 cm -3 and initial redshifts 2 = 10 (upper curves) and 2 = 20 
(lower curves). The horizontal dashed lines indicate the CMB 
temperature at these redshifts. Four different treatments of H2 
cooling are compared: our full model, CF1 (solid line); two vari- 
ants of this model, one which omits H+ and e _ collisions (CF2; 
dashed line), and one which omits H+, e — and He collisions (CF3; 
dot-dashed line); and the treatment used in most pr evious studies, 
which includes only collisions with H, but uses thelGalli fc Pallal 
il998l ) cooling rate rather than the IWrathmall fc Flowed d2007f ) 
rate. 



3.2 Choice of H2 cooling function 

To explore the sensitivity of the thermal evolution of pri- 
mordial gas, and in particular of the ability of the gas to 
cool to temperatures at which HD cooling dominates, to 
uncertainties and omissions in the treatment of H2 cooling, 
we ran a number of free-fall collapse and isobaric evolution 
models using different H2 cooling functions. Our reference 
model (hereafter CF1) uses the cooling function outlined in 
i|2.3l and used elsewhere in this paper; it includes the effects 
of collisions with H, H2, He, H + and e _ , and an ortho-para 
ratio that was computed self-consistently with the evolu- 
tion of the gas. We also examined the effects of omitting 
the H + and e~ contributions (CF2), and of omitting H + , 
e~ and He (CF3); note that in the latter case, the H2 cool- 
ing function essentially consists only of the H2-H contribu- 
tion, as the H2 fractions in our calculations are too small 
for H2-H2 collisions to ever become important. Finally, we 
examine the e ffect of includi n g only the H2-H contribution, 
but using the [G alli & Palla (1998) cooling rate instead of 
the lWrathmall fc Flower] (|2007l ) rate (CF4). Note that case 
CF4 assumes a fixed ortho-para ratio of 3:1, while the other 
treatments determine the ortho-para ratio self-consistently, 
as outlined above. However, the results of the previous sec- 
tion demonstrate that in practice this should not be a major 
source of error. 

In Figure [6] we show the temperature evolution of the 
gas as a function of density in two free-fall collapse models 
with an initial density rii — 1 cm -3 and redshifts z = 10 



and z = 20. The solid, dashed, dot-dashed and dotted lines 
correspond to CF1, CF2, CF3 and CF4 respectively. There 
are two important points to note about this plot. 
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Figure 7. Evolution of H2 and HD abundances as a function of 
gas number density in free-fall collapse models with initial den- 
sity ni = 1 cm -3 and initial redshift z = 20. We compare the 
behaviour for four different treatments of the H2 cooling: CF1 
(solid line), CF2 (dashed line), CF3 (dash-dotted line) and CF4 
(dotted line). Full details of these treatments are given in the text. 



First, it is clear that the temperature evolution of the 
gas in case CF1 differs from that in the other models over 
the whole range of density studied here. Although the ini- 
tial cooling of the gas is rapid in this model, this only lasts 
until the temperature reaches T ~ 300 K. At lower tem- 
peratures, the cooling of the gas slows down dramatically, 
allowing cooling in the other models to catch up and surpass 
it. The gas reaches a minimum temperature of 81 K in the 
z = 20 run and 53 K in the z = 10 run, significantly higher 
than the limits set by the CMB. 

Second, the temperature evolution of the gas in mod- 
els CF2, CF3 and CF4 differs noticeably at densities n < 
10 4 cm -3 and temperatures T < 500 K. The behaviour of 
models CF2 and CF4 is surprisingly similar, given the dif- 
ference in physical content of these two models, but the be- 
haviour of model CF3 is clearly different, with the gas in the 
latter model cooling less rapidly than in the other two. At 
densities higher than 10 4 cm -3 , however, all three models 
converge, and by the end of the simulation, the tempera- 
tures differ by no more than 10%. We obtain very similar 
results for free-fall collapse models with m = 0.03 cm -3 and 
m = 30 cm" 3 . 

The rapid cooling of the gas at early times in case CF1 is 
an obvious consequence of our inclusion of the cooling arising 
from H2-H + and H2-e~ collisions, but the relatively slow rate 
of cooling at later times (compared to the other models) 
at first seems somewhat counterintuitive: by adding extra 
coolants, we have made the gas cool more slowly! However, 
this puzzle is easy to solve if we examine the evolution of the 
H2 fraction in these simulations. In Figure UJ we show how 
the H2 and HD abundances evolve for models CF1-CF4; for 
clarity, we plot only the z = 20 case, although the behaviour 
in the z = 10 case is very similar. 

Figure [7] shows that significantly less H2 is produced 
in case CF1 than in the other cases, with the difference 
amounting to more than a factor of two at late times. This is 



a direct result of the rapid cooling of the gas at early times 
in run CF1. Most of the H2 that forms in all of the runs does 
so at early times, while the fractional ionization of the gas is 
still large. Enhanced cooling during this period reduces the 
rate of H~ formation (owing to the positive temperature de- 
pendence of the reaction coefficient for the formation of H~ 
by radiative association, reaction 1), and also increases the 
destruction rate of H~ by mutual neutralization (since the 
range increases for decreasing T, regardless of which partic- 
ular rate coefficient we adopt). The net result is a reduction 
in the H2 formation rate during this critical early period, 
and hence a reduction in the H2 abundance at late times. 

Figure [JJ also demonstrates that this reduction in the 
H2 abundance leads to a corresponding reduction in the HD 
abundance, which is a simple consequence of the fact that 
in most circumstances, reactions 39 and 41 dominate the 
production and destruction of HD, and are in equilibrium, 
implying that ihd oc ih 2 • 

The differences between runs CF2, CF3 and CF4 are 
also easy to u nderstand. At low d ensities and low tem- 
peratures, the iGalli fc Pallal l|l998h cooling function pro- 
vides significantly m o re co oling per H2 molecule than the 
IWrathmall fc Flowed (120071 ') . as Figure [2] demonstrates, and 
so gas in run CF4 can more easily reach a low tempera- 
ture than gas in run CF3. Below a temperature of about 
150 K, however, HD cooling begins to dominate in both 
models. As the H2 abundance does not differ greatly be- 
tween the two models, HD cooling is comparably effective 
in each, and the thermal evolution of the gas becomes insen- 
sitive to the choice of H2 cooling function. The temperature 
curves therefore begin to converge, with this convergence be- 
coming complete by the time that the gas reaches a density 
n ~ 10 4 cm -3 . At very high densities (n > 10 7 cm" 3 ), the 
gas once again becomes too warm for HD cooling to domi- 
nate, as the enhancement of the HD abundance by chemical 
fractionation becomes much less pronounced. However, at 
these densities, H2 is in LTE, and the only uncertainty in 
the H2 cooling rate comes from the small uncertainties in 
the energies of the various rotational and vibrational levels, 
and in the radiative transition rates. 

Including helium, as in run CF2, increases the H2 
cooling rate, particularly at low temperatures, and so 
the gas cools faster . How ever, the combination of the 
IWrathmall fc Flowed l|2007h rate for H2-H cooling with the 
H2-He cooling rate presented in >i2.3.3l scaled by the ap- 
propriate He:H ratio, coincide ntally results i n a to tal cool- 
ing rate that is similar to the IGalli fc Pallal j 19981 ) cooling 
rate: the two differ by no more than 25% in the temperature 
range 210 < T < 1000 K (assuming an ortho-para ratio of 
3: 1), despite the large disp arity in the H2-H coo l ing ra tes 
of IGalli fc Pallal (I1998T ) and IWrathmall fc Flowed (|2007f ) at 
these temperatures. 

The isobaric evolution models tell a similar story. Fig- 
ure [S^, shows how the gas temperature evolves with time in 
two representative models with initial density m = 1 cm -3 
and redshifts z — 10 (lower curves) and z = 20 (upper 
curves). We again find that the choice of H2 cooling func- 
tion affects the temperature evolution. Models CF2 and CF4 
again barely differ, while model CF3 differs from them sig- 
nificantly only at low temperatures (T < 700 K). At high 
temperatures (T S> 700K), cooling in all three of these mod- 
els is dominated the vibrational excitation of H2 by hydro- 
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gen atoms, and the vibra tional rates differ little betw een the 
iGalli fc Pallal (|l998f > and lWrathmall fc Flower! i|2007l ) treat- 
ments. On the other hand, if we include the effects of H + 
and e~ collisional excitation of H2, as in model CF1, we see 
a more substantial difference in the temperature evolution 
of the gas, at all temperatures T < 7000 K. Because of the 
high initial ionization, H2-H + collisions dominate, and the 
gas cools much faster than in the other models. We find sim- 
ilar results in our rii — 0.03 cm -3 and n\ — 30 cm -3 models. 

The importance of the difference in cooling time is un- 
clear, but likely depends upon the other relevant timescales 
in the problem. For instance, in the present example, the 
time required for the gas to cool to Tomb is much less than 
the Hubble time tn regardless of which treatment of H2 cool- 
ing is used. On the other hand, in our ri\ — 0.03cm -3 , z = 20 
model, illustrated in Figure [8]d, we find that t coo \ ~ tn for 
CF1, CF2 and CF4, but is approximately twice as long as 
the Hubble time if we use treatment CF3. 

We close this section by noting that despite the dif- 
ferences in the temperature evolution brought about by a 
change in H2 cooling function, in every case we have ex- 
amined the gas remains able to cool below the temperature 
reachable by H2 cooling alone. In other words, HD cooling 
is important in every case considered here, and the gas al- 
ways reaches the regime in which HD cooling dominates. 
However, in the free-fall models, the minimum temperature 
reached by the gas is always higher than the temperature 
floor set by the CMB, and varies depending on whether or 
not we include the effects of H2-H + and and H2-e~ collisions 
when computing the H2 cooling ra te. If, as has been hypoth - 
esized by some authors (see e.g. I Johnson fc Bromml l2006). 
gravitational fragmentation of the gas occurs only once the 
gas reaches its minimum temperature, then simulations that 
do not include these processes will produce fragments that 
are too small by roughly a factor of two. On the other hand, 
if the outcome of the fragmentation process is determined in 
part by the gas dynamics at earlier times, then the inclusion 
of these processes could conceivably enhance fragmentation, 
owing to the reduction in the cooling time of the gas at 
early times, and the much flatter temperature dependence 
of the H2-H + cooling rate compared to the H2-H cooling 
rate. Our very simple dynamical models do not allow us to 
explore these issues in any greater detail, but will hopefully 
motivate further work on the subject. 

3.3 Uncertainties in the reaction rate coefficients 



In £12.11 we discussed the large uncertainties that exist in 
some of the rate coefficients for reactions included in our 
chemical model. The most uncertain rates in our model are 
the destruction of H~ by associative detachment with H 
(reaction 2; see t ]2.1.1[) . the destruction of H~ by mutual 
neutralization with H + (reaction 5; again see i]2.1.1[) . the 
destruction of H2 by charge transfer with H + (reaction 7; 
see £!2.1.2p . and the three-body formation of H2 (reactions 
30 and 31; see i|2.1.5fl . In the following sections, we examine 
the individual effects of each of these uncertainties, before 
concluding by placing limits on the combined effect of all 
four uncertainties. 

Although we have examined the effects of these uncer- 
tainties for all of the combinations of n\ and z considered 
previously in this paper, for simplicity (and for clarity in 
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Figure 8. (a) As Figure[6] but showing the temperature evolution 
versus time for a set of models in which the evolution of the gas 
is isobaric. The initial density was n; = 1 cm -3 , and results are 
plotted for both 2 = 10 (lower set of lines) and 2 = 20 (upper set 
of lines). Four different treatments of H2 cooling are compared: 
CF1 (solid line), CF2 (dashed line), CF3 (dash-dotted line) and 
CF4 (dotted line). The horizontal dashed lines indicate the CMB 
temperature at redshifts z. = 10 and z = 20. For reference, the 
Hubble time at 2 = 10 is t H ~ 2.3 X 10 16 s and at 2 = 20 is 
*H — 9-0 X 10 15 s, where we have adopted th e standard WMAP3 
cosmological parameters llSpereel et alj|2007l) . (b) As (a), but for 
isobaric models with rij = 0.03 cm -3 ; note that only the 2 = 20 
case is plotted. 



the figures) we restrict our discussion here to one particular 
case: ni = 1 cm -3 and z = 20. Unless otherwise noted, we 
find very similar results for all of the other combinations of 
redshift and density that we have studied. 



3.3.1 Associative detachment and mutual neutralization 

In gas cooling and recombining from an initially ionized 
state, the amount of H2 formed is sensitive to the ratio be- 
tween the destruction rate of H~ by associative detachment 
with H (reaction 2) and by mutual neutralization with H + 
(reaction 5). An increase in k-z or a decrease in fes leads 
to associative detachment becoming the dominant destruc- 
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tion process at earlier times, when the fractional ioniza- 
tion of the gas is larger, and hence leads to a larger fi- 
nal H2 fraction; conversely, a decrease in or an increase 
in fcs means that mutual neutralization dominates for a 
longer period, and hence t he fin al H2 fraction is smaller. As 
iGlover. Savin fc Jappsenl (|2006t ) have already shown, in the 
absence of a substantial ultraviolet background, the effect 
on the final H2 fraction is not as large as might be feared: 
an order of magnitude change in ft2 or fcs alters the final 
H2 fraction by no more than a factor of a few. Nevertheless, 
this is enough to alter the temperature evolution of the gas 
by an appreciable amount, as we can see from Figures [5] and 

ESI 

Figure [U^ shows the effect on a representative free- 
fall collapse model of varying the associative detachment 
rate while keeping the mutual neutralization rate fixed. 
We plot result s from a model using our default value for 
ki, taken from ISchmeltekopf et all If 9671 ). and from mod- 
els using 'maximal' and 'minim al' values for ki taken from 
IGlover. Savin fc Jappsenl (|2006t ). The uncertainty in k^ in- 
troduces an uncertainty into the temperature evolution that 
persists throughout the simulation. Gas in simulations with 
a high value for (and hence higher H2 fractions) has a 
systematically lower temperature than the gas in simula- 
tions with a low value for ki. The difference between the 
simulations is particularly pronounced for densities in the 
range f 6 < n< 10 7 cm . The relatively rapid increase in 
the gas temperature at these densities occurs because HD 
reaches local thermodynamic equilibrium, and hence can no 
longer cool the gas so effectively. The causes the gas to be- 
gin heating up, which in turn reduces the HD abundance 
(as fractionation becomes less effective), causing the gas to 
warm further. The rate at which this process occurs de- 
pends upon the initial HD abundance, and hence on the H2 
abundance; reheating occurs more slowly when the H2 abun- 
dance is large. The gas temperature in this density regime 
can therefore differ by a factor of two or more, depending 
on which value is chosen for &2- 

In Figure [9p, we show how the same variation in k-z af- 
fects the temperature evolution in a representative isobaric 
model. In this case, the reduction in the H2 fraction resulting 
from a decrease in systematically delays cooling relative 
to our reference calculation. Similarly, an increase in £?2 ac- 
celerates cooling. In the present context, what is perhaps 
most interesting is the time taken to reach the temperature 
floor set by the CMB. This occurs after t~4x 1 14 s in the 
model with the largest value of fe2 and after t ~ 6 x 1 14 s 
in the model with the smallest value of ki- We find a sim- 
ilar degree of uncertainty in the cooling times in our other 
isobaric models. 

In Figure If 01 we examine the effect of varying the mu- 
tual neutralization rate while keeping the associative de- 
tachment rate fixed. We plot results fro m models per f ormed 
using mutual neutralization rates from Croft et al. I jl999l ) 
- our default - a s well as from iMoselev et al.l ( 1 97C ) and 
iDalgarno fc Leppl l|f987h . Figure EQ] shows that varying the 
mutual neutralization rate has very similar effects to vary- 
ing the associative detachment rate, except that the sense 
of the effect is reversed: a decrease in fcs has a similar effect 
to an increase in ki and vice versa. The size of the uncer- 
tainty introduced into the temperature evolution of the gas 
is comparable at low densities, and somewhat larger at high 



h o 
3 ° 
« 2 

ft 

B 



- 1 1 1 

: 

: 


1 1 1 1 ■ 1 1 1 1 1 1 rTTTT] 1 rTTTTTTj 1 1 1 lllll| 1 1 1 1 llll| 1 1 1 1 IIH 

: 

; 








' ' /' 




' /' 























10 100 1000 10 4 10 5 10 6 
Number density n (cm -3 ) 



10' 10° 



"T" 




10 1 



10 1 



10 l 



10' 



Time (s) 



Figure 9. (a) Temperature evolution as a function of gas den- 
sity for a free-fall collapse model with nj = 1 cm~ 3 and z = 
20. Three different values are used for the H - associative de- 
tachment rate coefficient (&2): our default value, taken from 
ISchmeltekopf et al.l ( |l96 7) (solid line); and 'maximal' and 'min- 
imal' values (dashed and dash-dotte d lines, respectively) taken 
from lGlover. Savin fc Japp sen (2006). The horizontal dashed line 
indicates the CMB temperature at z = 20. (b) As (a), but show- 
ing the time evolution of temperature in an isobaric model with 
the same initial conditions. 



densities, where the gas temperatures differ by as much as 
a factor of four. 

Nevertheless, although both rate coefficient uncertain- 
ties clearly affect the cooling of the gas, in neither case do 
they substantially change the outcome of the simulations. 
The gas still cools to temperatures low enough for chemical 
fractionation to significantly enhance HD, and so in each 
case HD cooling becomes dominant, further cooling the gas. 
In our free-fall collapse models, the minimum temperature 
reached by the gas does depend on the values of ki and fc 5 , 
and this may affect the characteristic fragment mass scale, 
although we would still expect any fragments to be smaller 
than in the case where only H2 cooling is effective. In our iso- 
baric models, the same minimum temperature is reached in 
every case, but the time taken to arrive there differs by up to 
a factor of two. The mass accretion dependent collapse times 
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Figure 10. (a) Temperature evolution as a function of gas density 
for a free-fall collapse model with n; = 1cm -3 and z = 20. Three 
different values are used for the H — + H + mutual neutralization 
rate coeffici ent (fc s ): ou r default rate (solid line; CDG99), taken 
from ICroft et al.l dl999|), along with a lar ge rate (dashed line; 
MOS70), taken from iMoselev et alj lll970h. and a smaller rate 
(dash-dotted line; DL87), taken from Dalgarno fc Leppl Jl987t) ■ 



The horizontal dashed line indicates the CMB temperature at 
2 = 20. (b) As (a), but showing the time evolution of temperature 
in an isobaric model with the same initial conditions. 
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Figure 11. (a) Temperature evolution as a function of gas density 
for a free-fall collapse model with n; = lcm~ 3 and z = 20. Three 
different values are used for the rate coefficient for H2 destruction 
by charge transfer with H + (reaction 7) . The rate used in our ref- 
erence model (solid line) is taken from Savin et al. (2004), but wc 
also show the effects of using rates from IShapiro fc Kand (119871 ) 
(dashed line) and lAbel et al l dl997r) (dash-dotted line). The hor- 
izontal dashed line indicates the CMB temperature at z = 20. 
(b) As (a), but showing the time evolution of temperature in an 
isobaric model with the same initial conditions. 



found bv lYoshida et all (|2003l ) suggest that such a factor of 
two uncertainty can be relevant, but further investigation 
requires a proper three-dimensional hydrodynamical treat- 
ment; our highly simplified dynamical models can take us 
no further. 

3.3.2 H2 charge transfer 

In Figure 1111 we examine the impact of varying the rate 
coefficient for H2 destruction by charge transfer (kr) in the 
context of representative free-fall collapse and isobaric evo- 
lution models. In the free-fall model, the effect of increasing 
k-j is to enable the gas to cool to lower temperatures. At 
first sight, this seems counterintuitive: by destroying H2, we 
make the gas colder. However, the key is that charge transfer 
is only an effective destruction mechanism at high temper- 



atures. By increasing the charge transfer rate, we delay the 
onset of rapid H2 cooling, and so when the gas does become 
able to cool rapidly, the fractional ionization is lower, and 
the cooling time is longer. Consequently, the gas remains 
warm for a longer period, forms more H~ (owing to the tem- 
perature dependence of reaction 1), and hence forms more 
H2. Nevertheless, the effect is relatively small: the uncer- 
tainty in the H2 fraction once the gas has cooled is no more 
than 5%, and the uncertainty in the minimum temperature 
is no more than 10%. 

In the isobaric model, we again see that the effect of in- 
creasing kj is to delay the onset of efficient H2 cooling. How- 
ever, once the gas begins cooling, the temperature evolution 
becomes convergent and the final outcome of the simulations 
is insensitive to the value of k-j. 
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3.3.3 Three-body H 2 formation 

Since we have assumed, following iPalla. Salpeter. fc Stahlerl 
(|l983h . that the rate coefficients for reactions 30 and 31 are 
related by £31 = &30/8, we can explore the effects of the 
uncertainty in the three-body rates simply by varying £30. 
In Figure \TZ\ we examine the effect of the uncertainty in re- 
action 30 in two representative models: one free-fall collapse 
model and one isobaric model. 

In the free-fall collapse model, the effect of the uncer- 
tainty is apparent only for densities n > 5 x I0 6 cm -3 . Be- 
tween this density and n = 10 8 cm -3 , the gas temperature 
increases slightly faster in the simulation that uses the larger 
iFlower fc H arris (2007) r ate than in t h e simu lation that uses 
the smaller Abel. Bryan, fc Normanl ([2002) rate, owing to 
the greater three-body H2 formation heating rate in the for- 
mer c ase. At n > 10 8 cm ~ 3 , however, the greater heating rate 
in the lFlower fc Harris! run is more than counterbalanced by 
the greater cooling provided by the larger abundance of H2, 
and s o the gas temperature inc reases at a slower rate than 
in the I Abel. Bryan, fc Normanl run. This difference persists 
until we terminate the simulation at n — 10 12 cm -3 , and the 
final temperatures differ by about 65%. 

In the isobaric model, there are no obvious differences 
between the two simulations. This is to be expected, as the 
gas density in these simulations never exceeds a few thou- 
sand particles per cubic centimeter, and so the three-body 
H2 formation rate remains extremely small throughout both 
simulations. 

In neither case does the uncertainty affect the cooling of 
the gas; in the free-fall model, it merely affects how quickly 
the gas subsequently re-heats. Therefore, whatever its im- 
pact on later stages of the star formation process, from the 
point of view of understanding the role and effectiveness of 
HD cooling, this particular source of uncertainty is clearly 
unimportant. 



3.3.4 Combining the uncertainties 

We close this discussion by examining two limiting cases, 
illustrated in Figure 1131 where we have selected the val- 
ues of the various uncertain rate coefficients in order to 
maximize and to minimize the degree of cooling. In our 
'maximal' model, we used a value of 5.0 x 10~ 9 cm 3 s _1 
for the H~ associative de tachment rate coefficient 
|Glover. Savin fc Jap pscn 200(]), and used values for the H~ 
mutual neutralization, H2 charge transfer and H2 thr e e-bod y 
formation r ate co effici ents from iDalgarno fc Leppl (1 1987). 
lAbel et all (|l997l ) and I Abel. Bryan, fc Normanl l|2002l ) re- 
spectively. In our 'minimal' model, we used a value of 
0.65 x 10~ 9 cm 3 s" 1 for the H~ associative detachmen t 
rate coefficient (again from[Gi over. Savin fc Jappscn 2006), 
along with values for the ot h er three rate coeffi c ients taken 
from iMoselev et all (ll970D . IShapiro fc Kand (|l987l ) and 
IFlower fc Harris (12007 ) respe c tively . Note that in both mod- 
els. the lWrathmall fc Flower! (|2007h H2 cooling rate is used 
and the effects of cooling from H2-H2, H2-He, H2-H + and 
H2-e~ collisions arc included; i.e. we do not couple the chem- 
ical rate uncertainties to the cooling rate uncertainties. 

In Figure If 3a . we examine the effects of the 'mini- 
mal' and 'maximal' models in the context of a represen- 
tative free-fall collapse calculation with an initial redshift 
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Figure 12. (a) Temperature evolution as a function of gas den- 
sity in two simulations with n; = 30 cm -3 and z = 20. Two 
different values are used for the rate coefficient for three-body 
H2 formation (£30): our defau lt value (solid line), taken from 
lAbel. Bryan, fc Normanl [|2002h, and a much larger value (dashed 
line) recently computed bv IFlower fc H arris (20q3). The horizon- 
tal dashed line indicates the CMB temperature at z = 20. (b) As 
(a), but showing the time evolution of temperature in an isobaric 
model with the same initial conditions. Note that in this plot, the 
results of the two simulations are indistinguishable. 



2 = 20. The temperature evolution in the two models dif- 
fers significantly at all densities, particularly in the interval 
10 4 < n < I0 8 cm . The minimum temperature reached 
by the gas in the 'maximal' model is 69.8 K, and in the 
'minimal' model is 141.5 K, a factor of two difference (cor- 
responding to a difference in the Jeans mass at this point 
of almost a factor of three). It is clear that in the 'mini- 
mal', HD cooling is of limited importance: although some 
HD cooling does occur, it never becomes completely domi- 
nant, and the gas temperature remains well above the floor 
set by the CMB. If the rates adopted in this model do prove 
to be the most accurate ones, then this would imply that HD 
cooling in collapsing gas may be significantly less effective 
than previously thought. 

Finally, in Figure 113b . we examine the effects of the 
'minimal' and 'maximal' models in the context of an isobaric 
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Figure 13. (a) Temperature evolution as a function of gas density 
in two simulations with n; = 1 cm -3 and z = 20, with values for 
the uncertain rate coefficients selected so as to minimize (solid 
line) or maximize (dashed line) the amount of cooling from H2. 
The horizontal dashed line indicates the CMB temperature at 
2 = 20. (b) As (a), but showing the time evolution of temperature 
in an isobaric model with the same initial conditions. 



evolution model with an initial redshift z = 20. Again, we see 
that the rate coefficient uncertainties significantly affect the 
cooling time of the gas, lengthening it in this case by a factor 
of three. However, they do not change the final outcome of 
the simulation: the gas still cools down to T ~ Tomb, and 
even for the 'minimal' model, the time required to reach 
T"cmb when m = 1 cm -3 is significantly less than a Hubble 
time. 



4 SUMMARY 

We have detailed the effect of several possible sources of un- 
certainty on the thermal evolution of primordial gas cooling 
from an initially hot and ionized state. These considered po- 
tential sources of uncertainty are the sensitivity of the low 
temperature H2 cooling rate to the ratio of ortho-H2 to para- 
H2 (which may differ significantly at low temperatures from 
the 3:1 ratio that is usually assumed), the continuing uncer- 



tainty in the form of the low temperature cooling rate for 
H2 excited by collisions with atomic hydrogen, the neglect 
by most previous authors of the contributions made to the 
H2 cooling rate by collisions of H2 with He, H + and e~ , and 
the large uncertainties that exist in the rates of several of 
the reactions responsible for determining the H2 fraction in 
the gas. 

We find that the first of these sources of uncertainty 
is unimportant. The standard assumption of a 3:1 ortho- 
para ratio is reasonably accurate at temperatures where H2 
cooling is effective. It only becomes significantly inaccurate 
for temperatures T < 100 K, but at these temperatures H2 
cooling is unimportant in comparison to HD cooling, and 
so this inaccuracy has no effect on the thermal evolution of 
the gas. However, the presented rates and formalism should 
prove useful when detailed complete emission models from 
cooling primordial gas are to be constructed. 

The sec ond source of uncert ainty - i.e. whether one uses 
the popular iGalli fc Pallal (|l998l ) parameterization to repre- 
sent th e cooling rate due to H2-H collisions, or the newer 
rates of IWrathmall fc Flower] (|2007l ) - can give the impres- 
sion of being important if one assumes that this process 
dominates the to tal H2 cooling rate. At low temperatures 
(T < 1000K), the lWrathmall fc Flower! <|2007l ) cooling func- 
tion p rovides significantly less cooling than the lGalli fc Pallal 
(1998) cooling function, and the use of the former in place 
of the latter lengthens the cooling time of the gas by about 
a factor of two. 

However, this simple comparison overstates the effect 
of this uncertainty. The problem stems from the assumption 
that H2-H collisions domina te. In practice, if one adopts the 
IWrathmall fc Flowerl \200ti ) cooling rate, then at low tem- 
peratures, H2-He collisions are more important, despite the 
low abundance of He relative to H. If one includes the effects 
of these collisions in the calculations, then the sensitivity of 
the outcome to the uncertainty in the H2-H rate becomes 
considerably less. 

A more important source of error in previous investiga- 
tions of HD formation in relic H 11 regions and other similar 
environments is the neglect of the effects of cooling due to 
collisions between H2 and protons and electrons. At early 
times during the cooling of the gas, the fractional ionization 
remains large, and H2-H + collisions dominate. When their 
effects are included, the cooling time of the gas is signifi- 
cantly decreased. However, a side effect of this rapid cooling 
is that less H2 forms, owing to the temperature dependence 
of the H2 formation rate. Because of this, less HD is formed, 
and so HD cooling is less effective. Therefore, the gas re- 
mains significantly warmer at later times than it would in 
models that did not include cooling from H2-H + collisions, 
although it remains cooler than it would be if HD cooling 
were not included. The increase in the minimum temper- 
ature in our free-fall models corresponds to an increase in 
the minimum Jeans mass of a factor of two, suggesting that 
previous studies of population III star formation in formerly- 
ionized regions may have underestimated the characteristic 
mass of the s t ars th at form. Moreover, previous work by 
lYoshida et alj l|2003h has shown that variations in the cool- 
ing time of the gas at low densities can have a pronounced ef- 
fect on the supply of cold gas available for population III star 
formation, and on t he tim ing of the collapse. Furthermore, 
lO'Shea fc Normanl (|2007h have shown that the amount of 
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H2 formed at low densities, which is sensitive to the ther- 
mal evolution of the gas, can affect the eventual accretion 
rate of gas onto the protostar, and hence may also affect its 
final mass. We therefore anticipate that the inclusion of the 
effects of H2-H + cooling may lead to clear differences in the 
outcome of such three-dimensional studies. 

Our investigation into the effects of the chemical rate 
coefficient uncertainties has shown that the large uncertain- 
ties in the associative detachment and mutual neutralization 
rates also have a significant impact on the thermal evolution 
of the gas. Although HD becomes the dominant coolant in 
all of the models considered, variations in the associative 
detachment and mutual neutralization rates alter the mini- 
mum temperature reached by the gas in our free-fall collapse 
simulations, and have a particularly pronounced effect on 
the temperature evolution at the end of the period during 
which HD cooling dominates. It is therefore quite plausi- 
ble that these uncertainties may also modify the outcome of 
multi-dimensional hydrodynamical simulations. 

On the other hand, our results display only a small sen- 
sitivity to the uncertainty in the H2 charge transfer reaction, 
other than at very early times in the isobaric runs, due to 
the fact that most of the H2 that forms in the gas does so 
at temperatures where this reaction is ineffective. Further 
investigation of the effects of this uncertainty may be of in- 
terest, but is clearly not a high priority. 

Finally, our results demonstrate that the evolution of 
gas in the HD-cooled regime is insensitive to the large un- 
certainty that exists in the three-body H2 formation rate, 
as three-body processes are unimportant at the densities at 
which HD dominates. However, this uncertainty does be- 
come important at densities n > 10 8 cm -3 . Its effect on the 
hydrodynamics of the gas at these densities remains uncer- 
tain. 

In summary, we discussed a number of hitherto ne- 
glected physical processes. The most crucial finding is that 
the molecular hydrogen cooling initiated by collisions with 
electrons, protons and neutral helium cannot be neglected 
in general, and should be included in future studies of the 
dynamics of primordial gas. 
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APPENDIX A: CHEMICAL NETWORK 

In Table [AT] we list the chemical reactions included in our 
model of primordial gas, along with the rate coefficients 
adopted and the references from which these rates were 
taken. For further details on some of the reactions, see also 

am 
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Table Al. List of reactions included in our chemical model 



No. Reaction Rate coefficient (cm 3 s x ) Ref. 

1 H + c-^H-+7 fci = dex[— 17.845 + 0.762 log T T ^ 6000 K 1 

+ 0.1523(log T) 2 

- 0.03274(logT) 3 ] 

= dcx[-16.4199 + 0.1998(logT) 2 T > 6000 K 

- 5.447 x 10~ 3 (logT) 4 
+ 4.0415 x 10- 5 (logT) 6 ] 

2 H~ + H -> H 2 + e~ See text 

3 H + H+^H++7 fc 3 = dex[-19.38 - 1.523 log T 2 

+ 1.118(logT) 2 - 0.1269(logT) 3 ] 

4 H + H+^H 2 +H+ fe 4 = 6.4 x 10- 10 3 

5 H" + H+ ->■ H + H See text 

6 H++e-^H + H fc 6 = 1.0xl0- 8 T ^ 617 K 4 

= 1.32 x io- 6 T-°' 76 T > 617 K 

7 H 2 +H+->H+ + H k 7 = [-3.3232183 x 10~ 7 5 

+ 3.3735382 x 10- 7 log T 

- 1.4491368 x 10" 7 (logT) 2 
+ 3.4172805 x 10- 8 (logT) 3 
-4.7813720 x 10" 9 (logT) 4 
+ 3.9731542 x 10~ 10 (logT) 5 

- 1.8171411 x 10 _11 (logT) 6 
+ 3.5311932 x 10- 13 (logT) 7 ] 
xexp( - 212 T 3715 ) 

8 H 2 +c- -+H + H + e- k 8 = 4.49 x 10- 9 T°- n cxp (-±2«58 ~j v = 6 

= 1.91 x 10- 9 T 0136 cxp (-£2 4 £Li) LTE 6 

9 H 2 +H-+H + H + H fe 9 = 6.67 x 10- 12 T°- 5 cxp [-(1 + v=0 7 

= 3.52 X 10- 9 cxp (-iSm) LTE 8 

10 H 2 +H 2 ^H 2 + H + H fcio = ( i.^re^~M-C)Td88i exp(-£™) V = Q g 

= 1.3 x 10~ 9 exp (-S222&) LTE 10 

11 H 2 +Hc^H + H + Hc fen =dcx [-27.029 + 3.801 log T- 22|87] 



dex 



v = 11 

-2.729 - 1.75 log T- 22474] LTE 11 



T J 

12 H + c- -+ H+ +c" +c" k 12 = cxp[-3.271396786 x 10 1 12 

+ 1.35365560 x 10 1 lnT e 

- 5.73932875 x 10°(lnT e ) 2 
+ 1.56315498 x 10°(lnT e ) 3 

- 2.87705600 x lO-^lnTe) 4 
+ 3.48255977 x 10" 2 (lnTe) 5 

- 2.63197617 x 10- 3 (lnT e ) 6 
+ 1.11954395 x 10- 4 (lnT e ) 7 

- 2.03914985 x 10- 6 (lnT e ) 8 ] 

13 H+ + e"- + H + 7 fe 13 , A = 1-269 xl0- 13 (2if 14 ) 1 - 503 Case A 13 

x [i. +(2Q 4 p)°- 470 ] -i.923 

fcis,B = 2.753 xl0- 14 (21fl4) 1 - 500 CaseB 13 

X [1.0+ (ilW88) - 407 ] -2.242 

14 H" +e- -> H + e" + e" k 14 = cxp[-l. 801849334 x 10 1 12 

+ 2.36085220 x 10°lnT e 

- 2.82744300 x 10- 1 (lnT e ) 2 
+ 1.62331664 x 10 _2 (lnT c ) 3 

- 3.36501203 x 10- 2 (lnT e ) 4 
+ 1.17832978 x 10" 2 (lnTe) 5 

- 1.65619470 x 10- 3 (lnT e ) 6 
+ 1.06827520 x 10" 4 (lnTe) 7 

- 2.63128581 x 10- 6 (lnT e ) 8 ] 
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Table Al - continued 



No. Reaction 



Rate coefficient (cm 3 s 1 ) 



Rcf. 



15 



19 



20 



21 
22 
23 
24 
25 
26 



H -> H + H + e" 



16 H+ + H-^H++c- 

17 He + e" -> He+ + e~ + e" 



18 He+ + ( 



He++ 



Hc+ 



He + 7 



He++ + e - — > He+ + 7 



H- 

H" 



+ H+ 
+ H+ 
H 2 + e" - 



^H 2 +H 
^H+H+H 
■ H~ +H 



H 2 + He+ -> He + H + H+ 



H 2 + He+ 



Hc^ 



H^ 



H -> He - 



+ He 
H+ 



fciB = 2.5634 x io-9T e 178186 
= exp[-2.0372609 x 10 1 
+ 1.13944933 x 10° lnT c 

- 1.4210135 x 10 _1 (lnT e ) 2 
+ 8.4644554 x 10- 3 (lnT o ) 3 

- 1.4327641 x 10" 3 (lnT o ) 4 
+ 2.0122503 x 10- 4 (lnT o ) 5 
+ 8.6639632 x 10- 5 (lnT e ) 6 
-2.5850097 x 10- 5 (lnT B ) 7 
+ 2.4555012 x 10- 6 (lnT c ) 8 

- 8.0683825 x 10- 8 (lnT c ) 9 ] 
: 6.9 x io- 9 T-°- 35 



fcl6 

kn 



fcl9,rr 
^19, rr 
fcl9,di 
kiO,A 
^20,B 

kn = 
kii = 
kiz = 

kl4 
k25 
k'26 



: 9.6 x io-'T" u au 
: exp[-4.409864886 x 10 1 
+ 2.391596563 x 10 1 lnT e 

- 1.07532302 x lO^lnTe) 2 
+ 3.05803875 x 10°(lnT e ) 3 

- 5.68511890 x 10- 1 (lnT e ) 4 
+ 6.79539123 x 10- 2 (lnT e ) s 

- 5.00905610 x 10- 3 (lnT e ) 6 
+ 2.06723616 x 10- 4 (lnT e ) 7 
-3.64916141 x 10- 6 (lnT e ) 8 ] 

: cxp[-6. 87104099 x 10 1 
+ 4.393347633 x lOMnTe 

- 1.84806699 x lO^lnTe) 2 
+ 4.70162649 x 10°(lnT c ) 3 

- 7.6924663 x 10-i(lnT e ) 4 
+ 8.113042 x 10~ 2 (lnTe) 5 

- 5.32402063 x 10- 3 (lnT e ) 6 
+ 1.97570531 x 10- 4 (lnT e ) 7 

- 3.16558106 x 10- 6 (lnT e ) 8 
A = io-HT" 5 [12.72 - 1.615 log T 

- 0.3162(log T) 2 + 0.0493(log T) 3 1 



= io-^T- 5 [11.19 - 1.676 logT 

- 0.2852(log T) 2 + 0.04433(log T) 3 
:1.9xlO- 3 T- 1 - 5 cxp(- 47 f2i) 
x [l.0 + 0.3cxp 

2.538 XlO" 13 (1262456 y 
2418500 \ 0.470 



( 94684 ^ 



X 



[1.0+ ( 



5.506 x 10 



T 
-14 I 



-0.5 



1262456 
T 

1.0+(«ffiM) - 407 

= 1.4X10"^ (^)-°- 
= 1.4 x 10- 7 (jgj) 
:2.7x 10- 8 T- 127 exp( 
:3.7x 10- 14 cxp(f ) 
: 7.2 x 10~ 15 
= 1.2 x 10-15 (^) 2B 



1.500 

-2.242 



43000 > 



Te^O.lcV 12 
T e > 0.1 cV 



T <C 8000 K 
T > 8000 K 



Case A 



Case B 



14 



12 



12 



Case A 15 
Case B 15 
Diclcctronic 16 



13 



13 



17 
17 
18 
19 
19 
20 
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Table Al - continued 



27 He + H+ ->■ Hc+ + H 

28 Hc++H-^He + H 

29 He + H" -> He + H + c" 

30 H + H + H^H 2 +H 

31 H + H + H 2 -» H 2 + H 2 

32 H + H + He ->■ H 2 + He 

33 D+ + c- -+D + 7 

34 D + H+ -> H + D+ 



35 H + D+ -> D + H 4 

36 H + D -> HD + 7 



37 H 2 + D -► HD + H 



38 
39 
40 

41 
42 

43 
44 
45 
46 
47 

48 
49 
50 
51 
52 

53 
54 
55 
56 
57 
58 
59 
60 
61 
62 



HD+ + H -> HD + H+ 
H 2 + D+ ->■ HD + H+ 
HD + H -> H 2 + D 

HD + H+ -» H 2 + D+ 
D + H+ — > HD+ + 7 
H + D+ — > HD+ + 7 
HD+ + e- ->• H + D 
D + c~ -> D+ + e~ +e 
He+ + D -> D+ + He 
He + D+ — > D + Hc+ 



H 2 + - 
HD H 



D 



HD+ + 
D + c- 
H + D- 
D + H" 
D + H" 
H + D- 
D + D" 
HD + c- 
HD + c- 
D 2 +c- 
H+ + D 
D+ + H 
D+ + D 



. HD+ +H 
— > HD + D+ 
^H++D 
D" +7 

► D + H~ 

► H + D- 

► HD + e" 

► HD + e" 

► D 2 + e" 
-^H + D- 
-^D + H- 
+ D + D" 

-» HD+ + c~ 
-> HD+ +c- 



fe 27 = 1.26 x 10- 9 T-°- 75 cxp ( 

= 4.0 x i - 37 T 4 - 74 
fc28 = 2.32 xl0- 7 (^)-°' 52 
fc 2g = 4.1 X 10- 17 T 2 cxp ( 
See text 
See text 

k 32 = 6.9 x 10-32T-0.4 

&33 = ki3 

fc 3 4 = 2.0xl0- 10 T o - 402 cxp(- 
- 3.31 x io- 17 T 148 



127500 
T , 



CXP I 22400 ) 
19870 \ 



T < 10000 K 
T > 10000 K 



^35 

&36 



&37 



&38 
&39 
&40 

fetl 
^42 
^43 

&45 
&47 



fc 4 g 

&52 
&53 

&55 
^56 
&57 
^58 
^59 
&60 

fcei 



: 3.44 x io- 10 T - 35 
: 2.06 xl0- lo T - 396 cxp (-f) 
+ 2.03 x io- 9 T-°' 332 
: 10~ 25 [2.80202 - 6.63697 InT 
+ 4.75619(lnT) 2 - 1.39325(lnT) 3 
+ 0.178259(ln T) 4 - 0.00817097(ln T) 5 ] 
: 10~ 25 exp[507.207 - 370.889 InT 
+ 104.854(lnT) 2 - 14.4192(ln T) 3 
+ 0.971469(ln T) 4 - 0.0258076(ln T) 5 ] 
: dcx [-56.4737 + 5.88886 log T 
+ 7.19692(log T) 2 
+ 2.25069(logT) 3 
- 2.16903(log T) 4 
+ 0.317887(log T) 5 ] 
:3.17x 10- 10 cxp(-5f 7 ) 

: [0.417 + 0.846 log T - 0.137(log T) 2 ] x 10" 



5.25 x 10- 



cxp 



4430 ' 
T , 
4430 
T 



173900 



5.25 x 10~ n exp 

1.1 x io-» exp (-4§8) 
: 3.9 x IO-1 9 18 exp (f) 

3-9 x 10-19 (&) 18 exp (*>) 

7.2 x 10- 8 T-°- 5 
k\i 

:1 .lxl0--(^)°- 25 



T < 2000 K 



T > 2000 K 



T ^ 200 K 
T > 200 K 



1.85 x 10- 9 T-°' 75 exp (- 



5.9 x 10 



-37^4.74 



127500 
T , 



T < 10000 K 

T > 10000 K 



(T_) 

\ 300 I 



0.002 



1.07 x lO- 9 
: 1.0 x lO" 9 exp (-^) 

6.4 x 10 — , 
6-4 x lO" 9 (^) - 
0.5 fc 2 
: 0.5 k 2 
k 2 

1.35 X 10" 9 T- 1 ' 27 cxp 
1.35 x 10" 9 T- 1 ' 27 cxp 
6.7 x 10- n T- 1:27 exp 

1.1 X 10-9 (^)-°- 4 
1.1 X 10-9 (^) 

1 . 3 xl0- 9 (^) 



cxp ( 



/_T_\0-41 
V 300 ) 



43000 

T 
43000 



-0.4 
-0.4 



21 

22 
23 



24 
25 



T < 2 x 10 5 K 26 
T > 2 x 10 5 K 

26 

10 < T < 200 K 27 
T > 200 K 



28 



25 
29 
30 

29 
31 

31 
32 
25 
31 
31 

33 
25 
34 
25 
34 
34 
35 
35 
25 
36 
36 
36 
31 
31 
31 



Uncertainties in H 2 and HD Chemistry and Cooling 



Table Al - continued 



63 
64 
65 
66 
67 
68 

69 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 



97 
98 



D" +e~ - 
D~ +H - 
D~ + He - 
D+ + H- 
H+ + D" 
D+ + D" 
H+ + D- 
H++D- 
HD+ + H 
HD+ + H 
HD+ + D 
HD+ + D 

d; 



D-2 
^2 

Hc+ + 
D + D 
D + H 



H- 
H" 
D" 
D- 
-D- 



I 2 + " 
H+ + D- 
HD+ +H 
HD+ +D 
HD+ +D 
D + D+ - 



H 

Do 



D 



■* D + e~ + e" 
D + H + e" 
-> D + He + e- 
-> D + H 
-> D + H 
-+ D + D 

^H 2 +D 
^H+H+D 

■ -» HD + H 

■ -> D + H + H 

■ -> HD + D 

■ -» D + H + D 
-> D 2 + H 

-> D + D- 
-+ D 2 +D 
-> D + D - 

->■ Hc + D 

D++ 7 

H 2 + D+ 

HD + H+ 
-> H 2 +D+ 
-^D++H 
-> D 2 +H+ 

D 2 + D+ 

D 2 +H+ 

HD+ + D 

HD + D+ 



- 11 



■D 



H 2 +D+ 
H 2 +D+ 



H^ 



- D 



HD+ + H 



^63 
ke4 
^65 

kea 
km 

&68 
&69 
&70 
k 71 

k72 

kr.i 

k74 
&75 

kn 
k 77 

k7S 
k79 

kso 
ksi 

^82 
^83 
fc 84 
&85 

k&e, 

k87 



kgo 
kgi 



= fel4 

= 1.5 x io- 17 T 2 cxp(-!2f20) 
= fc 5 
= fc 5 
= fc 5 

= 1.7xl0- 7 (^)-°- 5 

= 1.7xl0-'(&)-°- 5 

= 1.5xl0- 7 (^)_°- 5 

:L5 xlO- 7 (^) °' 5 

= 1.9X10"' 

= 1.9xl0- 7 (^)_°- 5 

:L5 xlO- 7 (^) °' J 

= 1.5xlO-'(&)-°- B 
= 2.0xl0-?(_r.)-°- 5 

= 2.0xl0- 7 (^)-°' 5 
= 3-03xl0-(^)-°- 52 
: 1.9 x 10- 19 T^- 8 



cxp (f ) 



CX P ( 22400 ) 



: fc 4 

: 1.0 X 10" 9 

: 1.0 X 10~ 9 

: 1.0 X 10~ 9 

: 1.0 X 10~ 9 

: ki 

- k^ 

:1.0XlO- 9 CXp(-^P) 
: 1.0 X 10~ 9 

-k 7 

= [1.04 x 10- 9 + 9.52 xlO- 9 ^) 
_l. 8 l x 10- 9 ( Ti r_) 2 ] cxp (-21000) 



92 


HD + H+ 


-> HD+ + H 


kg2 


— ^7 


93 


HD + H+ 


H+ +D 


^93 


- 1.0 X 10- 9 cxp( 21 6°°) 


94 


HD + D+ 


-» HD+ + D 


fc 94 


= k 7 


95 


HD + D+ 


-> D 2 + H+ 


^95 


= 1.0 x 10~ 9 


96 


HD + D+ 


^D++H 


^96 


= [3.54 x 10~ 9 +7.50 x 10" 



D 2 + H+ 
D 2 + H+ 



HD + D+ 
HD+ + D 



&97 
&98 



V 1UUUU / 

_ 2 . 92x 10-10 (_T_) 2 j exp (_2nm) 

= 2.1 x 10- 9 cxp(-fi) 
= [5.18 x 10- 11 + 3.05 x 10- 9 (1^0) 
" 5-42 x IO-10 (^j cxp (-20100) 



25 
25 
31 
25 
25 
25 
31 
31 
31 
31 
31 
31 
31 
31 
31 
31 
31 
31 
25 
37 
37 
38 
37 
25 
25 
38 
37 
25 
39 

25 
37 
25 
38 
39 



38 
39 



99 


D 2 4 


-H+ -> D+ +H 


&99 = 


= k 7 




25 


100 


D 2 -I 


- D+ -> D+ +D 


fcioo 


= k 7 




25 


101 


HD - 


f He+ -» HD+ + He 


feioi 


= k 2 5 




25 


102 


HD - 


f He+ -> He + H+ + D 


&102 


= 1.85 x 10~ 14 cxp ( 


35 

T 1 


35 


103 


HD - 


f He+ -> He + H + D+ 


fel03 


= 1.85 x 10- 14 cxp ( 


35 \ 
T ) 


35 


104 


D 2 i 


- He+ -<• D+ + He 


fcl04 


= 2.5 x 10~ 14 




38 


105 


D 2 4 


- He+ -<• He + D+ + D 


fcl05 


= 1.1 x 10- 13 T,-°' 24 


38 


106 


HD + D -> D 2 + H 


fcl06 


= 1.15 x 10-n exp ( 


3220 s ! 
T ) 


30 



28 S. C. 0. Glover & T. Abel 

Table Al - continued 



107 



D 2 + H -» HD + D 



108 
109 
110 
111 

112 
113 
114 
115 



HD 
HD 
HD 
HD 

D 2 4 
D 2 4 
D 2 4 
D 2 4 



•H - 
H 2 
•He 



H - 
-H 2 
-He - 



>H4 
-> H 
-+ H 
-> H 

D + 

•» D - 
-* D H 
-+DH 



D + 
+ D 
+ D- 
+-D- 

D + 

hD^ 
- D H 
-DH 



fc 107 = dex [-86.1558 + 4.53978 log T 
+ 33.5707(log T) 2 
- 13.0449(logT) 3 
+ 1.22017(logT) 4 
+ 0.0482453(log T) 5 ] 
= 2.67 x HT 10 exp (-^) 
H See text 

- H 2 See text 

- He See text 

-e" k ln = 5.09 x 10" 9 T° 128 exp 
= 1.04 x 10- <J T 218 cxp 
H fcii2 = fcg 

H 2 fcll3 = fcio 

He fcn4 = kn 



T I 
53070.7 N 



fell5 



: 8.24 x 10" 
■■ 2.75 x 10" 



exp 



105388 \ 
T 

53339.7 



T < 2200 K 28 



T > 2200 K 



LTE 



LTE 



40 

25 
25 
25 
6 



Note: T and T are the gas temperature in units of K and eV respectively. References are to the 

primary source o f data for each r eaction. 

References: 1: IWishard (Il979h. 2: iRamaker fe Peek! lll976h . 3: iKarpas. Anicich fc Huntress! 
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